{"id":300,"date":"2020-12-14T10:50:23","date_gmt":"2020-12-14T10:50:23","guid":{"rendered":"http:\/\/laserphotonics.uk\/?p=300"},"modified":"2020-12-14T10:50:23","modified_gmt":"2020-12-14T10:50:23","slug":"accuracy-of-bessel-functions-in-matlab","status":"publish","type":"post","link":"https:\/\/laserphotonics.uk\/?p=300","title":{"rendered":"Accuracy of Bessel Functions in MATLAB"},"content":{"rendered":"\n<p>by\u00a0PAVEL HOLOBORODKO\u00a0on\u00a0<abbr title=\"2016-05-12\">MAY 12, 2016<\/abbr><\/p>\n\n\n\n<figure class=\"wp-block-embed-wordpress wp-block-embed is-type-wp-embed is-provider-multiprecision-computing-toolbox-for-matlab\"><div class=\"wp-block-embed__wrapper\">\n<blockquote class=\"wp-embedded-content\" data-secret=\"4KvOjxscZQ\"><a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/\">Accuracy of Bessel Functions in MATLAB<\/a><\/blockquote><iframe title=\"&#8220;Accuracy of Bessel Functions in MATLAB&#8221; &#8212; Multiprecision Computing Toolbox for MATLAB\" class=\"wp-embedded-content\" sandbox=\"allow-scripts\" security=\"restricted\" style=\"position: absolute; clip: rect(1px, 1px, 1px, 1px);\" src=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/embed\/#?secret=4KvOjxscZQ\" data-secret=\"4KvOjxscZQ\" width=\"600\" height=\"338\" frameborder=\"0\" marginwidth=\"0\" marginheight=\"0\" scrolling=\"no\"><\/iframe>\n<\/div><\/figure>\n\n\n\n<p>Contents&nbsp;[<a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#\">hide<\/a>]<\/p>\n\n\n\n<ul><li><a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#Introduction\">1&nbsp;Introduction<\/a><\/li><li><a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#MATLAB_2016a\">2&nbsp;MATLAB 2016a<\/a><\/li><li><a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#Open_Source_Libraries\">3&nbsp;Open Source Libraries<\/a><ul><li><a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#GNU_GSL_21\">3.1&nbsp;GNU GSL 2.1<\/a><\/li><li><a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#Boost_1600\">3.2&nbsp;Boost 1.60.0<\/a><\/li><li><a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#Cephes_28\">3.3&nbsp;Cephes 2.8<\/a><\/li><li><a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#Numeric_Recipes_2007\">3.4&nbsp;Numeric Recipes, 2007<\/a><\/li><\/ul><\/li><li><a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#Commercial_Libraries\">4&nbsp;Commercial Libraries<\/a><ul><li><a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#NAG_Toolbox_for_MATLAB\">4.1&nbsp;NAG Toolbox for MATLAB<\/a><\/li><li><a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#Microsoft_CC_Language_and_Standard_Libraries\">4.2&nbsp;Microsoft C\/C++ Language and Standard Libraries<\/a><\/li><li><a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#Intel_Fortran\">4.3&nbsp;Intel Fortran<\/a><\/li><\/ul><\/li><li><a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#Conclusion\">5&nbsp;Conclusion<\/a><\/li><li><a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#References\">6&nbsp;References<\/a><\/li><\/ul>\n\n\n\n<h2>Introduction<\/h2>\n\n\n\n<p>In previous posts we studied accuracy of computation of modified Bessel functions:&nbsp;<a href=\"http:\/\/www.advanpix.com\/2016\/01\/05\/rational-approximations-for-the-modified-bessel-function-of-the-second-kind-k1-for-computations-with-double-precision\/\" target=\"_blank\" rel=\"noreferrer noopener\"><code>K<sub>1<\/sub>(x)<\/code><\/a>,&nbsp;<a href=\"http:\/\/www.advanpix.com\/2015\/11\/25\/rational-approximations-for-the-modified-bessel-function-of-the-second-kind-k0-for-computations-with-double-precision\/\" target=\"_blank\" rel=\"noreferrer noopener\"><code>K<sub>0<\/sub>(x)<\/code><\/a>,&nbsp;<a href=\"http:\/\/www.advanpix.com\/2015\/11\/11\/rational-approximations-for-the-modified-bessel-function-of-the-first-kind-i0-computations-double-precision\/\" target=\"_blank\" rel=\"noreferrer noopener\"><code>I<sub>0<\/sub>(x)<\/code><\/a>&nbsp;and&nbsp;<a href=\"http:\/\/www.advanpix.com\/2015\/11\/12\/rational-approximations-for-the-modified-bessel-function-of-the-first-kind-i1-for-computations-with-double-precision\/\" target=\"_blank\" rel=\"noreferrer noopener\"><code>I<sub>1<\/sub>(x)<\/code><\/a>. In spite of the fact that modified Bessel functions are easy to compute (they are monotonous and do not cross x-axis) we saw that MATLAB provides accuracy much lower than expected for&nbsp;<code>double<\/code>&nbsp;precision. Please refer to the pages for more details.<\/p>\n\n\n\n<p>Today we will investigate how accurately MATLAB computes Bessel functions of the first and second kind&nbsp;<code>Y<sub>n<\/sub>(x)<\/code>&nbsp;and&nbsp;<code>J<sub>n<\/sub>(x)<\/code>&nbsp;in double precision. Along the way, we will also check accuracy of the commonly used open source libraries.<\/p>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/bessel_j.png\" alt=\"BesselJ Functions (1st Kind, n=0,1,2)\"\/><\/figure>\n\n\n\n<p>Having extended precision routines, accuracy check of any function needs only few commands:<\/p>\n\n\n\n<pre class=\"wp-block-preformatted\">% Check accuracy of sine using 1M random points on (0, 16]\n&gt;&gt; mp.Digits(34);\n&gt;&gt; x = 16*(rand(1000000,1));\n&nbsp;\n&gt;&gt; f = sin(x);       % call built-in function for double precision sine\n&gt;&gt; y = sin(mp(x));   % compute sine values in quadruple precision\n&gt;&gt; z = abs(1-f.\/y);  % element-wise relative error\n&nbsp;\n&gt;&gt; max(z.\/eps)       % max. relative error in terms of 'eps'\nans = \n    0.5682267303295349594044472141263213\n&nbsp;\n&gt;&gt; -log10(max(z))    % number of correct decimal digits in result\nans = \n    15.89903811472552788729580391380739<\/pre>\n\n\n\n<p>Computed sine values have ~<code>15.9<\/code>&nbsp;correct decimal places and differ from \u2018true\u2019 function values by ~<code>0.5eps<\/code>. This is the highest accuracy we can expect from double precision floating-point arithmetic.<br><\/p>\n\n\n\n<p>The recent revision of&nbsp;<code>IEEE 754<\/code>&nbsp;standard<sup>[<a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#iee754\">1<\/a>]<\/sup>&nbsp;recommends all elementary functions to have such highest accuracy and be correctly rounded: the result should be as if the calculation had been performed in infinite precision, and then rounded (see Table 9.1 of the standard).<em>&nbsp;Computing elementary functions with correct rounding is a fascinating topic on its own, please refer to pioneering works of Gal<sup>[<a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#gal1\">2<\/a>,&nbsp;<a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#gal2\">3<\/a>]<\/sup>, Ziv<sup>[<a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#ziv\">4<\/a>]<\/sup>, recent developments by Lef\u00e8vre, Muller, Zimmermann, et al.<sup>[<a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#worst\">5<\/a>\u2013<a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#crlibm\">8<\/a>]<\/sup>&nbsp;and excellent handbooks on the subject<sup>[<a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#markstein\">9<\/a>\u2013<a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#handbook\">11<\/a>]<\/sup>.<\/em><\/p>\n\n\n\n<p>Most modern programming languages follow the&nbsp;<code>IEEE 754-2008<\/code>&nbsp;recommendation and their standard libraries<sup>[<a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#cpp\">12<\/a>,<a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#c99\">13<\/a>,<a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#fortran\">14<\/a>]<\/sup>&nbsp;compute elementary functions with correct rounding and guaranteed accuracy as we see in test above.<\/p>\n\n\n\n<p>However special functions are not covered by the standard and software libraries for its computations are considered to be&nbsp;<em>faithfully-accurate<\/em>. Here we study what libraries in fact can be trusted in computing Bessel functions in double precision.<\/p>\n\n\n\n<h2>MATLAB 2016a<\/h2>\n\n\n\n<p>Basic accuracy check applied for Bessel function of the second kind \u2013&nbsp;<code>Y<sub>0<\/sub>(x)<\/code>:<\/p>\n\n\n\n<pre class=\"wp-block-preformatted\">% Check accuracy of bessely(0,x) using 1M random points on (0, 16]\n&gt;&gt; f = bessely(0,x);     % built-in double precision routine\n&gt;&gt; y = bessely(0,mp(x)); % 'true' values computed using quadruple precision\n&gt;&gt; z = abs(1-f.\/y);\n&nbsp;\n&gt;&gt; max(z.\/eps)           % max. relative error in terms of 'eps'\nans = \n    1366834.014340578752744356956505545\n&nbsp;\n&gt;&gt; -log10(max(z))        % number of correct decimal digits in result\nans = \n    9.517843996632803435430170368588769<\/pre>\n\n\n\n<p>The problematic areas can be located using error scatter plot. Red lines show the expected limits of relative error (for convenience, we added another interval to cover whole x axis):<\/p>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MATLAB-Y_0-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MATLAB-Y_0-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<p>The spikes of error clearly coincide with zeros of Bessel function (<code>0.8935, 3.9576, 7.0860, ...<\/code>&nbsp;). To test accuracy near zeros in more detail we use the following function:<\/p>\n\n\n\n<pre class=\"wp-block-preformatted\">function func_accuracy_check(func, x0, near_zero, n)\n%func_accuracy_check - Checks accuracy of 'func' near its zero or extremum.\n%       \n%  1. Routine searches for zero\/extremum of 'func' starting from 'x0'\n%     Extended precision is used to be sure zero\/extremum is accurate. \n%    \n%  2. Relative error of 'func' is evaluated using 'n' points in close \n%     vicinity of found zero\/extremum. \n%\n%  Parameters:\n%   func - function handle \n%   x0 - starting point for zero\/extremum search\n%   near_zero - true\/false, check accuracy near zero\/extremum\n%   n - number of test points.\n%\n%  Usage example. \n%    Find zero of Y0(x) near x = 2 and check accuracy of the MATLAB \n%    function bessely(0,x) near found zero:\n%\n%    mp.Digits(34);\n%    u = @(x) bessely(0,x);\n%    func_accuracy_check(u,mp(2),true,3)\n%\n%    Output:\n%    Accuracy check near zero, x = 8.9357696627916752e-01:\n%    x = 8.9357696627916761e-01\tcorrect digits =   1.5\trel. error = 1.40e+14eps\n%    x = 8.9357696627916772e-01\tcorrect digits =   0.5\trel. error = 1.50e+15eps\n%    x = 8.9357696627916783e-01\tcorrect digits =   1.0\trel. error = 4.13e+14eps\n&nbsp;\n% Find root\/extremum of 'func' using extended precision:\noptions = optimset('TolX',mp('eps')); \nif near_zero\n    r = fzero(func,mp(x0),options);\n    fprintf('Accuracy check near zero, x = %.16e:\\n',r);    \nelse\n    %g = @(x) -func(x);    \n    r = fminsearch(func,mp(x0),options);\n    fprintf('Accuracy check near extremum, x = %.16e:\\n',r);    \nend;\n&nbsp;\n% Correctly rounded zero\/extremum in double precision:\nx = double(r); \nfor k=1:n\n    % Get next floating-point value after x:\n    x = nextafter(x);\n&nbsp;\n    f = func(x);      % function value computed by MATLAB's built-in bessely    \n    y = func(mp(x));  % 'true' function value computed in extended precision\n&nbsp;\n    z = abs(1-f.\/y);  % relative accuracy\n    fprintf('x = %.16e\\tcorrect digits = %4.1f\\trel. error = %5.2geps\\n',x,-log10(z),z.\/eps);  \nend<\/pre>\n\n\n\n<p>The&nbsp;<code>nextafter(x)<\/code>&nbsp;function computes next floating-point value just after x<sup>[<a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#nextafter\">15<\/a>]<\/sup>&nbsp;towards infinity. We use it to generate test arguments in close neighborhood of given point.<\/p>\n\n\n\n<p>Accuracy of&nbsp;<code>bessely(0,x)<\/code>&nbsp;near first zero of&nbsp;<code>Y<sub>0<\/sub>(x)<\/code>:<\/p>\n\n\n\n<pre class=\"wp-block-preformatted\">&gt;&gt; mp.Digits(34);\n&gt;&gt; u = @(x) bessely(0,x);\n&gt;&gt; func_accuracy_check(u,mp(2),true,5)\nAccuracy check near zero, x = 8.9357696627916752e-01:\nx = 8.9357696627916761e-01  correct digits =  1.5  rel. error = 1.4e+14eps\nx = 8.9357696627916772e-01  correct digits =  0.5  rel. error = 1.5e+15eps\nx = 8.9357696627916783e-01  correct digits =  1.0  rel. error = 4.1e+14eps\nx = 8.9357696627916794e-01  correct digits =  2.1  rel. error = 3.5e+13eps\nx = 8.9357696627916805e-01  correct digits =  1.4  rel. error = 1.6e+14eps<\/pre>\n\n\n\n<p>Accuracy of&nbsp;<code>bessely(0,x)<\/code>&nbsp;near second extremum of&nbsp;<code>Y<sub>0<\/sub>(x)<\/code>:<\/p>\n\n\n\n<pre class=\"wp-block-preformatted\">&gt;&gt; func_accuracy_check(u,mp(3),false,5)\nAccuracy check near extremum, x = 5.4296810407941351e+00:\nx = 5.4296810407941356e+00  correct digits = 15.6  rel. error =  1.10eps\nx = 5.4296810407941365e+00  correct digits = 16.2  rel. error =  0.32eps\nx = 5.4296810407941374e+00  correct digits = 16.0  rel. error =  0.42eps\nx = 5.4296810407941383e+00  correct digits = 16.0  rel. error =  0.42eps\nx = 5.4296810407941392e+00  correct digits = 15.6  rel. error =  1.10eps<\/pre>\n\n\n\n<p>MATLAB delivers minimal error near extremum but completely losses accuracy near function zero.<\/p>\n\n\n\n<p>Same analysis can be easily repeated for&nbsp;<code>J<sub>0<\/sub>(x)<\/code>:<\/p>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MATLAB-J_0-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MATLAB-J_0-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<p>Accuracy of&nbsp;<code>besselj(0,x)<\/code>&nbsp;near first zero of&nbsp;<code>J<sub>0<\/sub>(x)<\/code>:<\/p>\n\n\n\n<pre class=\"wp-block-preformatted\">&gt;&gt; u = @(x) besselj(0,x);\n&gt;&gt; func_accuracy_check(u,mp(2),true,5)\nAccuracy check near zero, x = 2.4048255576957728e+00:\nx = 2.4048255576957733e+00  correct digits =  0.6  rel. error = 1.2e+15eps\nx = 2.4048255576957738e+00  correct digits =  1.0  rel. error = 4.7e+14eps\nx = 2.4048255576957742e+00  correct digits =  0.8  rel. error = 7.4e+14eps\nx = 2.4048255576957747e+00  correct digits =  1.9  rel. error = 5.9e+13eps\nx = 2.4048255576957751e+00  correct digits =  1.4  rel. error = 1.6e+14eps<\/pre>\n\n\n\n<p>Accuracy of&nbsp;<code>besselj(0,x)<\/code>&nbsp;near first extremum of&nbsp;<code>J<sub>0<\/sub>(x)<\/code>:<\/p>\n\n\n\n<pre class=\"wp-block-preformatted\">&gt;&gt; func_accuracy_check(u,mp(2),false,5)\nAccuracy check near extremum, x = 3.8317059702075123e+00:\nx = 3.8317059702075129e+00  correct digits = 15.8  rel. error =  0.71eps\nx = 3.8317059702075134e+00  correct digits = 15.8  rel. error =  0.71eps\nx = 3.8317059702075138e+00  correct digits = 15.5  rel. error =  1.30eps\nx = 3.8317059702075142e+00  correct digits = 15.8  rel. error =  0.71eps\nx = 3.8317059702075147e+00  correct digits = 15.9  rel. error =  0.53eps<\/pre>\n\n\n\n<p>Situation is the same \u2013&nbsp;<code>besselj(0,x)<\/code>&nbsp;computes&nbsp;<code>J<sub>0<\/sub>(x)<\/code>&nbsp;incorrectly near zeros.<\/p>\n\n\n\n<p><strong>In fact, all routines in MATLAB for computing Bessel functions of integer order suffer from accuracy loss in vicinity of function zeros. Error can be arbitrary high, even to the point when none of the digits are correct.<\/strong><\/p>\n\n\n\n<p>Below we just show relative error plots with obvious peaks near function zeros \u2013 indication of accuracy loss (detailed analysis for every function can be easily done using routine above):<\/p>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MATLAB-Y_1-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MATLAB-Y_1-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MATLAB-Y_2-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MATLAB-Y_2-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MATLAB-Y_3-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MATLAB-Y_3-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MATLAB-J_1-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MATLAB-J_1-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MATLAB-J_2-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MATLAB-J_2-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MATLAB-J_3-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MATLAB-J_3-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<p>Certainly this situation is absolutely unacceptable, and TMW should re-consider its library for computing Bessel functions.<\/p>\n\n\n\n<p><strong>More surprisingly, the issue is endemic for nearly all libraries for computing Bessel functions in double precision!<\/strong>&nbsp;Below we present error plots for some commonly used open source and commercial libraries. All of them (except one) suffer from the same accuracy loss near function zeros \u2013 error may become arbitrary high.<\/p>\n\n\n\n<h2>Open Source Libraries<\/h2>\n\n\n\n<h3>GNU GSL 2.1<\/h3>\n\n\n\n<p><em>GNU Scientific Library<sup>[<a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#gsl\">16<\/a>]<\/sup><\/em><\/p>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/GNU-GSL-Y_0-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/GNU-GSL-Y_0-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/GNU-GSL-Y_1-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/GNU-GSL-Y_1-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/GNU-GSL-Y_2-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/GNU-GSL-Y_2-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/GNU-GSL-Y_3-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/GNU-GSL-Y_3-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/GNU-GSL-J_0-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/GNU-GSL-J_0-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/GNU-GSL-J_1-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/GNU-GSL-J_1-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/GNU-GSL-J_2-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/GNU-GSL-J_2-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/GNU-GSL-J_3-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/GNU-GSL-J_3-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<h3>Boost 1.60.0<\/h3>\n\n\n\n<p><em>Free peer-reviewed portable C++ source libraries<sup>[<a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#boost\">17<\/a>]<\/sup><\/em><\/p>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/BOOST-Y_0-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/BOOST-Y_0-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/BOOST-Y_1-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/BOOST-Y_1-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/BOOST-Y_2-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/BOOST-Y_2-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/BOOST-Y_3-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/BOOST-Y_3-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/BOOST-J_0-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/BOOST-J_0-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/BOOST-J_1-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/BOOST-J_1-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/BOOST-J_2-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/BOOST-J_2-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/BOOST-J_3-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/BOOST-J_3-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<h3>Cephes 2.8<\/h3>\n\n\n\n<p><em>C and C++ language special functions math library<sup>[<a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#cephes\">18<\/a>]<\/sup><\/em><\/p>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/CEPHES-Y_0-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/CEPHES-Y_0-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/CEPHES-Y_1-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/CEPHES-Y_1-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/CEPHES-Y_2-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/CEPHES-Y_2-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/CEPHES-Y_3-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/CEPHES-Y_3-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/CEPHES-J_0-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/CEPHES-J_0-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/CEPHES-J_1-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/CEPHES-J_1-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/CEPHES-J_2-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/CEPHES-J_2-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/CEPHES-J_3-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/CEPHES-J_3-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<h3>Numeric Recipes, 2007<\/h3>\n\n\n\n<p><em>3rd Edition: The Art of Scientific Computing<sup>[<a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#press\">19<\/a>]<\/sup><\/em><\/p>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/Numeric-Recipes-Y_0-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/Numeric-Recipes-Y_0-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/Numeric-Recipes-Y_1-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/Numeric-Recipes-Y_1-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/Numeric-Recipes-Y_2-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/Numeric-Recipes-Y_2-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/Numeric-Recipes-Y_3-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/Numeric-Recipes-Y_3-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/Numeric-Recipes-J_0-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/Numeric-Recipes-J_0-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/Numeric-Recipes-J_1-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/Numeric-Recipes-J_1-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/Numeric-Recipes-J_2-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/Numeric-Recipes-J_2-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/Numeric-Recipes-J_3-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/Numeric-Recipes-J_3-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<h2>Commercial Libraries<\/h2>\n\n\n\n<h3>NAG Toolbox for MATLAB<\/h3>\n\n\n\n<p><em>Mark 24 \u2013 newest version<sup>[<a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#nag\">25<\/a>]<\/sup><\/em><\/p>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/NAG-Y_0-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/NAG-Y_0-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/NAG-Y_1-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/NAG-Y_1-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/NAG-J_0-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/NAG-J_0-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/NAG-J_1-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/NAG-J_1-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<h3>Microsoft C\/C++ Language and Standard Libraries<\/h3>\n\n\n\n<p><em>Part of several generations of MS Visual Studio&nbsp;<sup>[<a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#press\">21<\/a>]<\/sup><\/em><\/p>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MSVC-Y_0-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MSVC-Y_0-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MSVC-Y_1-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MSVC-Y_1-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MSVC-Y_2-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MSVC-Y_2-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MSVC-Y_3-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MSVC-Y_3-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MSVC-J_0-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MSVC-J_0-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MSVC-J_1-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MSVC-J_1-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MSVC-J_2-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MSVC-J_2-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MSVC-J_3-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/MSVC-J_3-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<h3>Intel Fortran<\/h3>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/INTEL-FORTRAN-Y_0-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/INTEL-FORTRAN-Y_0-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/INTEL-FORTRAN-Y_1-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/INTEL-FORTRAN-Y_1-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/INTEL-FORTRAN-Y_2-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/INTEL-FORTRAN-Y_2-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/INTEL-FORTRAN-Y_3-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/INTEL-FORTRAN-Y_3-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/INTEL-FORTRAN-J_0-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/INTEL-FORTRAN-J_0-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/INTEL-FORTRAN-J_1-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/INTEL-FORTRAN-J_1-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/INTEL-FORTRAN-J_2-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/INTEL-FORTRAN-J_2-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/INTEL-FORTRAN-J_3-I.png\" alt=\"\"\/><\/figure>\n\n\n\n<figure class=\"wp-block-image\"><img src=\"https:\/\/www.advanpix.com\/wp-content\/uploads\/2016\/05\/INTEL-FORTRAN-J_3-II.png\" alt=\"\"\/><\/figure>\n\n\n\n<h2>Conclusion<\/h2>\n\n\n\n<p>All tested commercial (<code>MATLAB, NAG, Intel<\/code>&nbsp;and&nbsp;<code>Microsoft<\/code>) and open source libraries (<code>GNU GSL, Boost, CEPHES<\/code>&nbsp;and&nbsp;<code>Numeric Recipes<\/code>) exhibit severe accuracy degradation near zeros of Bessel functions.&nbsp;<strong>In fact, computed function values have no correct digits in close vicinity of zeros.<\/strong><\/p>\n\n\n\n<p>Intel mathematical library delivers accurate results for&nbsp;<code>Y<sub>0<\/sub>(x), Y<sub>1<\/sub>(x), J<sub>0<\/sub>(x)<\/code>&nbsp;and&nbsp;<code>J<sub>1<\/sub>(x)<\/code>&nbsp;but suffers from the same accuracy loss for higher orders&nbsp;<code>n=2,3,...<\/code>.<\/p>\n\n\n\n<p>The difficulty of approximating Bessel functions near zeros is well known<sup>[<a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#capprox\">22<\/a>]<\/sup>&nbsp;and efficient algorithms were developed for both cases \u2013 double<sup>[<a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#harris\">23<\/a>]<\/sup>&nbsp;and arbitrary precision<sup>[<a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#mpfr\">24<\/a>]<\/sup>. We hope maintainers of the libraries will take this analysis into account and will provide fixes in future releases.<\/p>\n\n\n\n<p>Until then, please rely on high-accuracy tools, e.g. Advanpix toolbox which is the only way to compute Bessels correctly in&nbsp;<code>MATLAB<\/code>&nbsp;at the moment.<\/p>\n\n\n\n<h2>References<\/h2>\n\n\n\n<ol><li>IEEE Computer Society (2008),&nbsp;<a href=\"http:\/\/www.csee.umbc.edu\/~tsimo1\/CMSC455\/IEEE-754-2008.pdf\" target=\"_blank\" rel=\"noreferrer noopener\">IEEE Standard for Binary Floating-Point Arithmetic<\/a>.<\/li><li>S. Gal.&nbsp;<em>Computing elementary functions: A new approach for achieving high accuracy and good performance<\/em>. In Accurate Scientific Computations. Lecture Notes in Computer Science, volume 235, pages 1\u201316. Springer-Verlag, Berlin, 1986.<\/li><li>S. Gal and B. Bachelis.&nbsp;<em>An accurate elementary mathematical library for the IEEE floating point standard<\/em>. ACM Transactions on Mathematical Software, 17(1):26\u201345, March 1991.<\/li><li>A. Ziv,&nbsp;<em>Fast Evaluation of Elementary Mathematical Functions with Correctly Rounded Last Bit<\/em>, ACM Trans. Math. Software, vol. 17, no. 3, pp. 410-423, 1991.<\/li><li>V. Lefevre and J.-M. Muller,&nbsp;<em>Worst Cases for Correct Rounding of the Elementary Functions in Double Precision<\/em>, Proc. 15th IEEE Symp. Computer Arithmetic (ARITH 15), N. Burgess and L. Ciminiera, eds., pp. 111-118, 2001<\/li><li>D. Defour, G. Hanrot, V. Lefevre, J.-M. Muller, N. Revol, and P. Zimmermann,&nbsp;<em>Proposal for a Standardization of Mathematical Function Implementation in Floating-Point Arithmetic<\/em>, Numerical Algorithms,vol. 37, nos. 1-4, pp. 367-375, 2004<\/li><li>D. Stehle and P. Zimmermann.&nbsp;<em>Gal\u2019s accurate tables method revisited<\/em>. In Proceedings of the 17th IEEE Symposium on Computer Arithmetic (ARITH-17). IEEE Computer Society Press, Los Alamitos, CA, 2005.<\/li><li>Catherine Daramy-Loirat, David Defour, Florent de Dinechin, Matthieu Gallet, Nicolas Gast, and Jean-Michel Muller.&nbsp;<em>CR-LIBM: A library of correctly rounded elementary functions in double-precision<\/em>, February 29, 2009.<\/li><li>P. W. Markstein.&nbsp;<em>IA-64 and Elementary Functions: Speed and Precision<\/em>. Hewlett Packard Professional Books. Prentice Hall, 2000.<\/li><li>J.-M. Muller.&nbsp;<em>Elementary Functions, Algorithms and Implementation<\/em>. Birkh\u00e4user Boston, 2nd Edition, 2006<\/li><li>Muller, J.-M., Brisebarre, N., de Dinechin, F., Jeannerod, C.-P., Lef\u00e8vre, V., Melquiond, G., Revol, N., Stehl\u00e9, D., Torres, S..&nbsp;<em>Handbook of Floating-Point Arithmetic<\/em>. Birkh\u00e4user Boston, 2009.<\/li><li>ISO\/IEC 14882: Programming Language C++ (Latest revisions).<\/li><li>ISO\/IEC 9899: Programming language C (Latest revisions).<\/li><li>JTC1\/SC22\/WG5:&nbsp;<a href=\"http:\/\/www.nag.co.uk\/sc22wg5\/\" target=\"_blank\" rel=\"noreferrer noopener\">FORTRAN<\/a>.<\/li><li><a href=\"http:\/\/jp.mathworks.com\/matlabcentral\/newsreader\/view_thread\/192\" target=\"_blank\" rel=\"noreferrer noopener\">Next floating-point number in MATLAB.<\/a>&nbsp;See also C\/C++ standards.<\/li><li>M. Galassi et al,&nbsp;<a href=\"http:\/\/www.gnu.org\/software\/gsl\/\">GNU Scientific Library Reference Manual<\/a>, 3rd Edition (January 2009), ISBN 0954612078. Version 2.1, November 11, 2015.<\/li><li><a href=\"http:\/\/www.boost.org\/\" target=\"_blank\" rel=\"noreferrer noopener\">Boost \u2013 free peer-reviewed portable C++ source libraries<\/a>. Version 1.60.0, December 17th, 2015.<\/li><li><a href=\"http:\/\/www.moshier.net\/\" target=\"_blank\" rel=\"noreferrer noopener\">CEPHES C and C++ language special functions math library<\/a>. Version 2.8, November 4, 2014.<\/li><li>Press, William H. and Teukolsky, Saul A. and Vetterling, William T. and Flannery, Brian P. Numerical Recipes 3rd Edition: The Art of Scientific Computing, 2007, Cambridge University Press, New York, USA.<\/li><li><a href=\"https:\/\/www.gnu.org\/software\/libc\/\">GNU C Library<\/a>, February 19, 2016<\/li><li><a href=\"https:\/\/msdn.microsoft.com\/en-us\/library\/hh875057.aspx\">Microsoft C\/C++ Language and Standard Libraries<\/a>.<\/li><li>J. F. Hart, E. W. Cheney, C. L. Lawson, H. J. Maehly, C. K. Mesztenyi, J. R. Rice, H. G. Thatcher, and C. Witzgall.&nbsp;<em>Computer Approximations<\/em>. Robert E. Krieger, 1978.<\/li><li>J. Harrison.&nbsp;<em><a href=\"http:\/\/www.cl.cam.ac.uk\/~jrh13\/papers\/bessel.pdf\" target=\"_blank\" rel=\"noreferrer noopener\">Fast and Accurate Bessel Function Computation<\/a><\/em>, Proceedings of ARITH19, the 19th IEEE Conference on Computer Arithmetic, IEEE Computer Society Press, 2009, pp. 104-113.<\/li><li>Fousse, G. Hanrot, V. Lefevre, P. Pelissier, and P. Zimmermann, 2007,&nbsp;<a href=\"http:\/\/www.mpfr.org\/\" target=\"_blank\" rel=\"noreferrer noopener\"><em>MPFR: A multiple-precision binary floating-point library with correct rounding<\/em><\/a>, ACM TOMS, 33(2), 13:1\u201315<\/li><li>The NAG Toolbox for MATLAB, Mark 24 (the newest version).<\/li><\/ol>\n\n\n\n<p>{&nbsp;1&nbsp;comment\u2026 read it below or&nbsp;<a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#respond\">add one<\/a>&nbsp;}<strong>Stefan<\/strong><a href=\"https:\/\/www.advanpix.com\/2016\/05\/12\/accuracy-of-bessel-functions-in-matlab\/#comment-224062\">September 18, 2020 at 6:29 pm<\/a><\/p>\n\n\n\n<p>Absolutely amazing work!<br>This is a serious issue, that I run into with my students, when expanding functions in a Bessel series. I can\u2019t tell you how often I have said \u201cyes, that looks sh*t because Matlab\u2019s Bessel function are sh*t. We have to stop after a few terms in the series\u201d.<\/p>\n","protected":false},"excerpt":{"rendered":"<p>by\u00a0PAVEL HOLOBORODKO\u00a0on\u00a0MAY 12, 2016 Contents&nbsp;[hide] 1&nbsp;Introduction 2&nbsp;MATLAB 2016a 3&nbsp;Open Source Libraries 3.1&nbsp;GNU GSL 2.1 3.2&nbsp;Boost 1.60.0 3.3&nbsp;Cephes 2.8 3.4&nbsp;Numeric Recipes, 2007 4&nbsp;Commercial Libraries 4.1&nbsp;NAG Toolbox for MATLAB 4.2&nbsp;Microsoft C\/C++ Language and Standard Libraries 4.3&nbsp;Intel Fortran 5&nbsp;Conclusion 6&nbsp;References Introduction In previous posts we studied accuracy of computation of modified Bessel functions:&nbsp;K1(x),&nbsp;K0(x),&nbsp;I0(x)&nbsp;and&nbsp;I1(x). In spite of the fact [&hellip;]<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":[],"categories":[1],"tags":[],"_links":{"self":[{"href":"https:\/\/laserphotonics.uk\/index.php?rest_route=\/wp\/v2\/posts\/300"}],"collection":[{"href":"https:\/\/laserphotonics.uk\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/laserphotonics.uk\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/laserphotonics.uk\/index.php?rest_route=\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/laserphotonics.uk\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=300"}],"version-history":[{"count":1,"href":"https:\/\/laserphotonics.uk\/index.php?rest_route=\/wp\/v2\/posts\/300\/revisions"}],"predecessor-version":[{"id":301,"href":"https:\/\/laserphotonics.uk\/index.php?rest_route=\/wp\/v2\/posts\/300\/revisions\/301"}],"wp:attachment":[{"href":"https:\/\/laserphotonics.uk\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=300"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/laserphotonics.uk\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=300"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/laserphotonics.uk\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=300"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}