Last active
December 6, 2020 14:53
-
-
Save ntessore/0e1764986731589eb843ee1bda3e38ea to your computer and use it in GitHub Desktop.
compute the logarithm of Bessel functions of large argument and order
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
// compute the logarithm of Bessel functions of large argument and order | |
// | |
// author: Nicolas Tessore <n.tessore@ucl.ac.uk> | |
// date: 6 Dec 2020 | |
// | |
#include <math.h> | |
// compute the Bessel function ratio J_{nu-1}(x)/J_nu(x) from a continued | |
// fraction | |
// | |
// reference: Rothwell, Commun. Numer. Meth. Engng 2008; 24:237-249 | |
// | |
double rjnu(double nu, double x) | |
{ | |
int m; | |
double r, a, c, d, delt; | |
const double eps = 1e-15; | |
const double tiny = 1e-30; | |
r = 2*nu/x; | |
c = r; | |
d = 0; | |
for(m = 1;; ++m) | |
{ | |
nu += 1; | |
a = (1-2*(m&1))*2*nu/x; | |
c = a + 1/c; | |
if(c == 0) | |
c = tiny; | |
d = a + d; | |
if(d == 0) | |
d = tiny; | |
d = 1/d; | |
delt = c*d; | |
r = r*delt; | |
if(fabs(delt - 1) < eps) | |
break; | |
} | |
return r; | |
} | |
// compute Bessel function logarithm ln[J_{nu+n}(x)] from ln[J_nu(x)] | |
// | |
// reference: Rothwell, Commun. Numer. Meth. Engng 2008; 24:237-249 | |
// + explicit tracking of sign change if `sgn` is not NULL | |
// | |
double lnjnupn(double nu, int n, double x, double lnjnu, double* sgn) | |
{ | |
int k; | |
double r, s, c, y, t; | |
r = rjnu(nu+n, x); | |
k = 0; | |
// Kahan summation | |
s = lnjnu; | |
c = 0; | |
for(; n > 0; --n) | |
{ | |
if(r < 0) | |
k ^= 1; | |
y = -log(fabs(r)) - c; | |
t = s + y; | |
c = (t - s) - y; | |
s = t; | |
r = 2*(nu+n-1)/x - 1/r; | |
} | |
if(sgn) | |
*sgn = k ? -1 : +1; | |
return s; | |
} | |
// test program for spherical Bessel functions -- delete | |
#include <stdio.h> | |
int main() | |
{ | |
// argument | |
double x = 5.0; | |
// order of start value | |
double nu = 0; | |
// abs log of start value: Log[Abs[SphericalBesselJ[0, 5.0`20]]] | |
double lnjnu = -1.6513810824626862774; | |
// track sign change | |
double sgn; | |
// how many times to raise the order | |
int n = 10000; | |
// compute the result: note +1/2 order of spherical Bessel function | |
lnjnu = lnjnupn(nu+0.5, n, x, lnjnu, &sgn); | |
// logjnu now contains j(nu+n, z) and sgn is the sign change from -1 | |
printf("jnu = %+f exp(%.16f)\n", -sgn, lnjnu); | |
// the true value: Log[SphericalBesselJ[10000, 5.0`20]] | |
printf(" [%+f exp(%.16f)]\n", 1.0, -72950.74713290145977); | |
return 0; | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Note that the
lnjnupn
function is good to compute a single valueln[J_{nu+n}(x)]
fromln[J_nu(x)]
. If all valuesln[J_{nu+k}(x)]
fork = 1, ..., n
are desired, it is much more efficient to compute first all ratiosrk
fromr
and subsequently compute the partial sums forwards.