bessel_zero.c

多个常用的特殊函数的例子

/* specfunc/bessel_zero.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_airy.h>
#include <gsl/gsl_sf_pow_int.h>
#include <gsl/gsl_sf_bessel.h>

#include "error.h"

#include "bessel_olver.h"

/* For Chebyshev expansions of the roots as functions of nu,
 * see [G. Nemeth, Mathematical Approximation of Special Functions].
 * This gives the fits for all nu and s <= 10.
 * I made the fits for other values of s myself [GJ].
 */

/* Chebyshev expansion: j_{nu,1} = c_k T_k*(nu/2), nu <= 2 */
static const double coef_jnu1_a[] = {
  3.801775243633476,
  1.360704737511120,
 -0.030707710261106,
  0.004526823746202,
 -0.000808682832134,
  0.000159218792489,
 -0.000033225189761,
  0.000007205599763,
 -0.000001606110397,
  0.000000365439424,
 -0.000000084498039,
  0.000000019793815,
 -0.000000004687054,
  0.000000001120052,
 -0.000000000269767,
  0.000000000065420,
 -0.000000000015961,
  0.000000000003914,
 -0.000000000000965,
  0.000000000000239,
 -0.000000000000059,
  0.000000000000015,
 -0.000000000000004,
  0.000000000000001
};


/* Chebyshev expansion: j_{nu,1} = nu c_k T_k*((2/nu)^(2/3)), nu >= 2 */
static const double coef_jnu1_b[] = {
  1.735063412537096,
  0.784478100951978,
  0.048881473180370,
 -0.000578279783021,
 -0.000038984957864,
  0.000005758297879,
 -0.000000327583229,
 -0.000000003853878,
  0.000000002284653,
 -0.000000000153079,
 -0.000000000000895,
  0.000000000000283,
  0.000000000000043,
  0.000000000000010,
 -0.000000000000003
};


/* Chebyshev expansion: j_{nu,2} = c_k T_k*(nu/2), nu <= 2 */
static const double coef_jnu2_a[] = {
  6.992370244046161,
  1.446379282056534,
 -0.023458616207293,
  0.002172149448700,
 -0.000246262775620,
  0.000030990180959,
 -0.000004154183047,
  0.000000580766328,
 -0.000000083648175,
  0.000000012317355,
 -0.000000001844887,
  0.000000000280076,
 -0.000000000042986,
  0.000000000006658,
 -0.000000000001039,
  0.000000000000163,
 -0.000000000000026,
  0.000000000000004,
 -0.000000000000001
};


/* Chebyshev expansion: j_{nu,2} = nu c_k T_k*((2/nu)^(2/3)), nu >= 2 */
static const double coef_jnu2_b[] = {
  2.465611864263400,
  1.607952988471069,
  0.138758034431497,
 -0.003687791182054,
 -0.000051276007868,
  0.000045113570749,
 -0.000007579172152,
  0.000000736469208,
 -0.000000011118527,
 -0.000000011919884,
  0.000000002696788,
 -0.000000000314488,
  0.000000000008124,
  0.000000000005211,
 -0.000000000001292,
  0.000000000000158,
 -0.000000000000004,
 -0.000000000000003,
  0.000000000000001
};


/* Chebyshev expansion: j_{nu,3} = c_k T_k*(nu/3), nu <= 3 */
static const double coef_jnu3_a[] = {
  10.869647065239236,
   2.177524286141710,
  -0.034822817125293,
   0.003167249102413,
  -0.000353960349344,
   0.000044039086085,
  -0.000005851380981,
   0.000000812575483,
  -0.000000116463617,
   0.000000017091246,
  -0.000000002554376,
   0.000000000387335,
  -0.000000000059428,
   0.000000000009207,
  -0.000000000001438,
   0.000000000000226,
  -0.000000000000036,
   0.000000000000006,
  -0.000000000000001
};


/* Chebyshev expansion: j_{nu,3} = nu c_k T_k*((3/nu)^(2/3)), nu >= 3 */
static const double coef_jnu3_b[] = {
  2.522816775173244,
  1.673199424973720,
  0.146431617506314,
 -0.004049001763912,
 -0.000039517767244,
  0.000048781729288,
 -0.000008729705695,
  0.000000928737310,
 -0.000000028388244,
 -0.000000012927432,
  0.000000003441008,
 -0.000000000471695,
  0.000000000025590,
  0.000000000005502,
 -0.000000000001881,
  0.000000000000295,
 -0.000000000000020,
 -0.000000000000003,
  0.000000000000001
};


/* Chebyshev expansion: j_{nu,4} = c_k T_k*(nu/4), nu <= 4 */
static const double coef_jnu4_a[] = {
  14.750310252773009,
   2.908010932941708,
  -0.046093293420315,
   0.004147172321412,
  -0.000459092310473,
   0.000056646951906,
  -0.000007472351546,
   0.000001031210065,
  -0.000000147008137,
   0.000000021475218,
  -0.000000003197208,
   0.000000000483249,
  -0.000000000073946,
   0.000000000011431,
  -0.000000000001782,
   0.000000000000280,
  -0.000000000000044,
   0.000000000000007,
  -0.000000000000001
};


/* Chebyshev expansion: j_{nu,4} = nu c_k T_k*((4/nu)^(2/3)), nu >= 4 */
static const double coef_jnu4_b[] = {
  2.551681323117914,
  1.706177978336572,
  0.150357658406131,
 -0.004234001378590,
 -0.000033854229898,
  0.000050763551485,
 -0.000009337464057,
  0.000001029717834,
 -0.000000037474196,
 -0.000000013450153,
  0.000000003836180,
 -0.000000000557404,
  0.000000000035748,
  0.000000000005487,
 -0.000000000002187,
  0.000000000000374,
 -0.000000000000031,
 -0.000000000000003,
  0.000000000000001
};



/* Chebyshev expansion: j_{nu,5} = c_k T_k*(nu/5), nu <= 5 */
static const double coef_jnu5_a[] = {
  18.632261081028211,
   3.638249012596966,
  -0.057329705998828,
   0.005121709126820,
  -0.000563325259487,
   0.000069100826174,
  -0.000009066603030,
   0.000001245181383,
  -0.000000176737282,
   0.000000025716695,
  -0.000000003815184,
   0.000000000574839,
  -0.000000000087715,
   0.000000000013526,
  -0.000000000002104,
   0.000000000000330,
  -0.000000000000052,
   0.000000000000008,
  -0.000000000000001
};


/* Chebyshev expansion: j_{nu,5} = nu c_k T_k*((5/nu)^(2/3)), nu >= 5 */
/* FIXME: There is something wrong with this fit, in about the
 * 9th or 10th decimal place.
 */
static const double coef_jnu5_b[] = {
  2.569079487591442,
  1.726073360882134,
  0.152740776809531,
 -0.004346449660148,
 -0.000030512461856,
  0.000052000821080,
 -0.000009713343981,
  0.000001091997863,
 -0.000000043061707,
 -0.000000013779413,
  0.000000004082870,
 -0.000000000611259,
  0.000000000042242,
  0.000000000005448,
 -0.000000000002377,
  0.000000000000424,
 -0.000000000000038,
 -0.000000000000002,
  0.000000000000002
};


/* Chebyshev expansion: j_{nu,6} = c_k T_k*(nu/6), nu <= 6 */
static const double coef_jnu6_a[] = {
  22.514836143374042,
   4.368367257557198,
  -0.068550155285562,
   0.006093776505822,
  -0.000667152784957,
   0.000081486022398,
  -0.000010649011647,
   0.000001457089679,
  -0.000000206105082,
   0.000000029894724,
  -0.000000004422012,
   0.000000000664471,
  -0.000000000101140,
   0.000000000015561,
  -0.000000000002416,
   0.000000000000378,
  -0.000000000000060,
   0.000000000000009,
  -0.000000000000002
};


/* Chebyshev expansion: j_{nu,6} = nu c_k T_k*((6/nu)^(2/3)), nu >= 6 */
static const double coef_jnu6_b[] = {
  2.580710285494837,
  1.739380728566154,
  0.154340696401691,
 -0.004422028860168,
 -0.000028305272624,
  0.000052845975269,
 -0.000009968794373,
  0.000001134252926,
 -0.000000046841241,
 -0.000000014007555,
  0.000000004251816,
 -0.000000000648213,
  0.000000000046728,
  0.000000000005414,
 -0.000000000002508,
  0.000000000000459,
 -0.000000000000043,
 -0.000000000000002,
  0.000000000000002
};


/* Chebyshev expansion: j_{nu,7} = c_k T_k*(nu/7), nu <= 7 */
static const double coef_jnu7_a[] = {
  26.397760539730869,
   5.098418721711790,
  -0.079761896398948,
   0.007064521280487,
  -0.000770766522482,
   0.000093835449636,
  -0.000012225308542,
   0.000001667939800,
  -0.000000235288157,
   0.000000034040347,
  -0.000000005023142,
   0.000000000753101,
  -0.000000000114389,
   0.000000000017564,
  -0.000000000002722,
   0.000000000000425,
  -0.000000000000067,
   0.000000000000011,
  -0.000000000000002
};


/* Chebyshev expansion: j_{nu,7} = nu c_k T_k*((7/nu)^(2/3)), nu >= 7 */
static const double coef_jnu7_b[] = {
  2.589033335856773,
  1.748907007612678,
  0.155488900387653,
 -0.004476317805688,
 -0.000026737952924,
  0.000053459680946,
 -0.000010153699240,
  0.000001164804272,
 -0.000000049566917,
 -0.000000014175403,
  0.000000004374840,
 -0.000000000675135,
  0.000000000050004,
  0.000000000005387,
 -0.000000000002603,
  0.000000000000485,
 -0.000000000000047,
 -0.000000000000002,
  0.000000000000002
};


/* Chebyshev expansion: j_{nu,8} = c_k T_k*(nu/8), nu <= 8 */
static const double coef_jnu8_a[] = {
  30.280900001606662,
   5.828429205461221,
  -0.090968381181069,
   0.008034479731033,
  -0.000874254899080,
   0.000106164151611,
  -0.000013798098749,
   0.000001878187386,
  -0.000000264366627,
   0.000000038167685,
  -0.000000005621060,
   0.000000000841165,
  -0.000000000127538,
   0.000000000019550,
  -0.000000000003025,
   0.000000000000472,
  -0.000000000000074,
   0.000000000000012,
  -0.000000000000002
};


/* Chebyshev expansion: j_{nu,8} = nu c_k T_k*((8/nu)^(2/3)), nu >= 8 */
static const double coef_jnu8_b[] = {
  2.595283877150078,
  1.756063044986928,
  0.156352972371030,
 -0.004517201896761,
 -0.000025567187878,
  0.000053925472558,
 -0.000010293734486,
  0.000001187923085,
 -0.000000051625122,
 -0.000000014304212,
  0.000000004468450,
 -0.000000000695620,
  0.000000000052500,
  0.000000000005367,
 -0.000000000002676,
  0.000000000000505,
 -0.000000000000050,
 -0.000000000000002,
  0.000000000000002
};


/* Chebyshev expansion: j_{nu,9} = c_k T_k*(nu/9), nu <= 9 */
static const double coef_jnu9_a[] = {
  34.164181213238386,
   6.558412747925228,
  -0.102171455365016,
   0.009003934361201,
  -0.000977663914535,
   0.000118479876579,
  -0.000015368714220,
   0.000002088064285,
  -0.000000293381154,
   0.000000042283900,
  -0.000000006217033,
   0.000000000928887,
  -0.000000000140627,
   0.000000000021526,
  -0.000000000003326,
   0.000000000000518,
  -0.000000000000081,
   0.000000000000013,
  -0.000000000000002
};


/* Chebyshev expansion: j_{nu,9} = nu c_k T_k*((9/nu)^(2/3)), nu >= 9 */
static const double coef_jnu9_b[] = {
  2.600150240905079,
  1.761635491694032,
  0.157026743724010,
 -0.004549100368716,
 -0.000024659248617,
  0.000054291035068,
 -0.000010403464334,
  0.000001206027524,
 -0.000000053234089,
 -0.000000014406241,
  0.000000004542078,
 -0.000000000711728,
  0.000000000054464,
  0.000000000005350,
 -0.000000000002733,
  0.000000000000521,
 -0.000000000000052,
 -0.000000000000002,
  0.000000000000002
};


/* Chebyshev expansion: j_{nu,10} = c_k T_k*(nu/10), nu <= 10 */
static const double coef_jnu10_a[] = {
  38.047560766184647,
   7.288377637926008,
  -0.113372193277897,
   0.009973047509098,
  -0.001081019701335,
   0.000130786983847,
  -0.000016937898538,
   0.000002297699179,
  -0.000000322354218,
   0.000000046392941,
  -0.000000006811759,
   0.000000001016395,
  -0.000000000153677,
   0.000000000023486,
  -0.000000000003616,
   0.000000000000561,
  -0.000000000000095,
   0.000000000000027,
  -0.000000000000013,
   0.000000000000005
};


/* Chebyshev expansion: j_{nu,10} = nu c_k T_k*((10/nu)^(2/3)), nu >= 10 */
static const double coef_jnu10_b[] = {
  2.604046346867949,
  1.766097596481182,
  0.157566834446511,
 -0.004574682244089,
 -0.000023934500688,
  0.000054585558231,
 -0.000010491765415,
  0.000001220589364,
 -0.000000054526331,
 -0.000000014489078,
  0.000000004601510,
 -0.000000000724727,
  0.000000000056049,
  0.000000000005337,
 -0.000000000002779,
  0.000000000000533,
 -0.000000000000054,
 -0.000000000000002,
  0.000000000000002
};


/* Chebyshev expansion: j_{nu,11} = c_k T_k*(nu/22), nu <= 22 */
static const double coef_jnu11_a[] = {
  49.5054081076848637,
  15.33692279367165101,
 -0.33677234163517130,
  0.04623235772920729,
 -0.00781084960665093,
  0.00147217395434708,
 -0.00029695043846867,
  0.00006273356860235,
 -0.00001370575125628,
  3.07171282012e-6,
 -7.0235041249e-7,
  1.6320559339e-7,
 -3.843117306e-8,
  9.15083800e-9,
 -2.19957642e-9,
  5.3301703e-10,
 -1.3007541e-10,
  3.193827e-11,
 -7.88605e-12,
  1.95918e-12,
 -4.9020e-13,
  1.2207e-13,
 -2.820e-14,
  5.25e-15,
 -1.88e-15,
  2.80e-15,
 -2.45e-15
};


/* Chebyshev expansion: j_{nu,12} = c_k T_k*(nu/24), nu <= 24 */
static const double coef_jnu12_a[] = {
  54.0787833216641519,
  16.7336367772863598,
 -0.36718411124537953,
  0.05035523375053820,
 -0.00849884978867533,
  0.00160027692813434,
 -0.00032248114889921,
  0.00006806354127199,
 -0.00001485665901339,
  3.32668783672e-6,
 -7.5998952729e-7,
  1.7644939709e-7,
 -4.151538210e-8,
  9.87722772e-9,
 -2.37230133e-9,
  5.7442875e-10,
 -1.4007767e-10,
  3.437166e-11,
 -8.48215e-12,
  2.10554e-12,
 -5.2623e-13,
  1.3189e-13,
 -3.175e-14,
  5.73e-15,
  5.6e-16,
 -8.7e-16,
 -6.5e-16
};


/* Chebyshev expansion: j_{nu,13} = c_k T_k*(nu/26), nu <= 26 */
static const double coef_jnu13_a[] = {
  58.6521941921708890,
  18.1303398137970284,
 -0.39759381380126650,
  0.05447765240465494,
 -0.00918674227679980,
  0.00172835361420579,
 -0.00034800528297612,
  0.00007339183835188,
 -0.00001600713368099,
  3.58154960392e-6,
 -8.1759873497e-7,
  1.8968523220e-7,
 -4.459745253e-8,
  1.060304419e-8,
 -2.54487624e-9,
  6.1580214e-10,
 -1.5006751e-10,
  3.679707e-11,
 -9.07159e-12,
  2.24713e-12,