fix for scaled besseli at negative non-integer orders
Brian Gough
bjg at gnu.org
Tue Oct 28 09:25:03 CDT 2008
Hello,
In the current tip there is a problem with the scaling factor for
besseli with negative non-integer order in liboctave/lo-specfun.cc.
Computing the ratio of scaled and unscaled besseli(alpha,2) should
give the same value for all orders but gives a different result for
orders like -3.5:
$ ./run-octave
GNU Octave, version 3.1.51+
....
octave:1> for alpha = [1,2,3,3.5,-1,-2,-3,-3.5] ; r=besseli(alpha,2,1)/besseli(alpha,2); disp(r); end
0.13534 # exp(-2) correct ratio
0.13534
0.13534
0.13534
0.13534
0.13534
0.13534
8.6239 # should also be 0.13534 = exp(-2)
The current code in zbesi() computes besseli(-alpha,z) as
zbesi(alpha,z) + (2/pi)*sin(pi*alpha)*zbesk(alpha,z)
This is correct in the unscaled case, but needs an extra term to
compensate for the different scaling factors of zbesi and zbesk in the
scaled case.
The scaling factor of zbesi() is si=exp(-abs(real(z))) but for zbesk it
is sk=exp(z). Hence the formula in the scaled case needs to be
zbesi(alpha,z) + (si/sk)*(2/pi)*sin(pi*alpha)*zbesk(alpha,z)
i.e. an extra factor of exp(-z-abs(real(z))). I've attached three
patches which
1) add all the different scaling factors for besselj,y,i,k,h to the
documentation
2) add tests for the bessel functions for combinations of
odd,even,fractional order and real and complex argument
3) correct the calculation of zbesi for negative order with scaling.
I've put the same correction into cbesi as well, although the single
precision calls to besseli don't currently seem to work (they always
return NaN for me).
--
Brian Gough
GNU Scientific Library -
http://www.gnu.org/software/gsl/
-------------- next part --------------
# HG changeset patch
# User Brian Gough <bjg at gnu.org>
# Date 1225198390 14400
# Node ID 71422651b143cd95dd6e757541885a0d769e5084
# Parent b4a1a54721917756500c0ec0b9252e69808d72db
add all scaling factors to bessel documentation
diff -r b4a1a5472191 -r 71422651b143 scripts/specfun/bessel.m
--- a/scripts/specfun/bessel.m Mon Oct 27 14:36:25 2008 -0400
+++ b/scripts/specfun/bessel.m Tue Oct 28 08:53:10 2008 -0400
@@ -26,21 +26,23 @@
##
## @table @code
## @item besselj
-## Bessel functions of the first kind.
+## Bessel functions of the first kind. If the argument @var{opt} is supplied,
+## the result is multiplied by @code{exp(-abs(imag(x)))}.
## @item bessely
-## Bessel functions of the second kind.
+## Bessel functions of the second kind. If the argument @var{opt} is supplied,
+## the result is multiplied by @code{exp(-abs(imag(x)))}.
## @item besseli
-## Modified Bessel functions of the first kind.
+## Modified Bessel functions of the first kind. If the argument @var{opt} is supplied,
+## the result is multiplied by @code{exp(-abs(real(x)))}.
## @item besselk
-## Modified Bessel functions of the second kind.
+## Modified Bessel functions of the second kind. If the argument @var{opt} is supplied,
+## the result is multiplied by @code{exp(x)}.
## @item besselh
## Compute Hankel functions of the first (@var{k} = 1) or second (@var{k}
-## = 2) kind.
+## = 2) kind. If the argument @var{opt} is supplied, the result is multiplied by
+## @code{exp (-I*@var{x})} for @var{k} = 1 or @code{exp (I*@var{x})} for
+## @var{k} = 2.
## @end table
-##
-## If the argument @var{opt} is supplied, the result is scaled by the
-## @code{exp (-I*@var{x})} for @var{k} = 1 or @code{exp (I*@var{x})} for
-## @var{k} = 2.
##
## If @var{alpha} is a scalar, the result is the same size as @var{x}.
## If @var{x} is a scalar, the result is the same size as @var{alpha}.
diff -r b4a1a5472191 -r 71422651b143 src/DLD-FUNCTIONS/besselj.cc
--- a/src/DLD-FUNCTIONS/besselj.cc Mon Oct 27 14:36:25 2008 -0400
+++ b/src/DLD-FUNCTIONS/besselj.cc Tue Oct 28 08:53:10 2008 -0400
@@ -388,21 +388,23 @@
\n\
@table @code\n\
@item besselj\n\
-Bessel functions of the first kind.\n\
+Bessel functions of the first kind. If the argument @var{opt} is supplied, \n\
+the result is multiplied by @code{exp(-abs(imag(x)))}.\n\
@item bessely\n\
-Bessel functions of the second kind.\n\
+Bessel functions of the second kind. If the argument @var{opt} is supplied,\n\
+the result is multiplied by @code{exp(-abs(imag(x)))}.\n\
@item besseli\n\
-Modified Bessel functions of the first kind.\n\
+Modified Bessel functions of the first kind. If the argument @var{opt} is supplied,\n\
+the result is multiplied by @code{exp(-abs(real(x)))}.\n\
@item besselk\n\
-Modified Bessel functions of the second kind.\n\
+Modified Bessel functions of the second kind. If the argument @var{opt} is supplied,\n\
+the result is multiplied by @code{exp(x)}.\n\
@item besselh\n\
Compute Hankel functions of the first (@var{k} = 1) or second (@var{k}\n\
- = 2) kind.\n\
+= 2) kind. If the argument @var{opt} is supplied, the result is multiplied by\n\
+ at code{exp (-I*@var{x})} for @var{k} = 1 or @code{exp (I*@var{x})} for\n\
+ at var{k} = 2.\n\
@end table\n\
-\n\
-If the argument @var{opt} is supplied, the result is scaled by the\n\
- at code{exp (-I*@var{x})} for @var{k} = 1 or @code{exp (I*@var{x})} for\n\
- @var{k} = 2.\n\
\n\
If @var{alpha} is a scalar, the result is the same size as @var{x}.\n\
If @var{x} is a scalar, the result is the same size as @var{alpha}.\n\
-------------- next part --------------
# HG changeset patch
# User Brian Gough <bjg at gnu.org>
# Date 1225198441 14400
# Node ID 16699ae4ef0a65eed6127d239ee87239f2dd80df
# Parent 71422651b143cd95dd6e757541885a0d769e5084
add bessel function tests
diff -r 71422651b143 -r 16699ae4ef0a src/ChangeLog
--- a/src/ChangeLog Tue Oct 28 08:53:10 2008 -0400
+++ b/src/ChangeLog Tue Oct 28 08:54:01 2008 -0400
@@ -1,3 +1,7 @@
+2008-10-28 Brian Gough <bjg at gnu.org>
+
+ * DLD-FUNCTIONS/besselj.cc: Added tests.
+
2008-10-23 John W. Eaton <jwe at octave.org>
* oct-hist.c (initialize_history): New arg, read_history_file)
diff -r 71422651b143 -r 16699ae4ef0a src/DLD-FUNCTIONS/besselj.cc
--- a/src/DLD-FUNCTIONS/besselj.cc Tue Oct 28 08:53:10 2008 -0400
+++ b/src/DLD-FUNCTIONS/besselj.cc Tue Oct 28 08:54:01 2008 -0400
@@ -632,6 +632,288 @@
}
/*
+%! # Test values computed with GP/PARI version 2.3.3
+%!
+%!shared alpha, x, jx, yx, ix, kx, nix
+%!
+%! # Bessel functions, even order, positive and negative x
+%! alpha = 2; x = 1.25;
+%! jx = 0.1710911312405234823613091417;
+%! yx = -1.193199310178553861283790424;
+%! ix = 0.2220184483766341752692212604;
+%! kx = 0.9410016167388185767085460540;
+%!
+%!assert(besselj(alpha,x), jx, 100*eps)
+%!assert(bessely(alpha,x), yx, 100*eps)
+%!assert(besseli(alpha,x), ix, 100*eps)
+%!assert(besselk(alpha,x), kx, 100*eps)
+%!assert(besselh(alpha,1,x), jx + I*yx, 100*eps)
+%!assert(besselh(alpha,2,x), jx - I*yx, 100*eps)
+%!
+%!assert(besselj(alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
+%!assert(bessely(alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
+%!assert(besseli(alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
+%!assert(besselk(alpha,x,1), kx*exp(x), 100*eps)
+%!assert(besselh(alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
+%!assert(besselh(alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
+%!
+%!assert(besselj(-alpha,x), jx, 100*eps)
+%!assert(bessely(-alpha,x), yx, 100*eps)
+%!assert(besseli(-alpha,x), ix, 100*eps)
+%!assert(besselk(-alpha,x), kx, 100*eps)
+%!assert(besselh(-alpha,1,x), jx + I*yx, 100*eps)
+%!assert(besselh(-alpha,2,x), jx - I*yx, 100*eps)
+%!
+%!assert(besselj(-alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
+%!assert(bessely(-alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
+%!assert(besseli(-alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
+%!assert(besselk(-alpha,x,1), kx*exp(x), 100*eps)
+%!assert(besselh(-alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
+%!assert(besselh(-alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
+%!
+%! x *= -1;
+%! yx = -1.193199310178553861283790424 + 0.3421822624810469647226182835*I;
+%! kx = 0.9410016167388185767085460540 - 0.6974915263814386815610060884*I;
+%!
+%!assert(besselj(alpha,x), jx, 100*eps)
+%!assert(bessely(alpha,x), yx, 100*eps)
+%!assert(besseli(alpha,x), ix, 100*eps)
+%!assert(besselk(alpha,x), kx, 100*eps)
+%!assert(besselh(alpha,1,x), jx + I*yx, 100*eps)
+%!assert(besselh(alpha,2,x), jx - I*yx, 100*eps)
+%!
+%!assert(besselj(alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
+%!assert(bessely(alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
+%!assert(besseli(alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
+%!assert(besselk(alpha,x,1), kx*exp(x), 100*eps)
+%!assert(besselh(alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
+%!assert(besselh(alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
+%!
+%! # Bessel functions, odd order, positive and negative x
+%! alpha = 3; x = 2.5;
+%! jx = 0.2166003910391135247666890035;
+%! yx = -0.7560554967536709968379029772;
+%! ix = 0.4743704087780355895548240179;
+%! kx = 0.2682271463934492027663765197;
+%!
+%!assert(besselj(alpha,x), jx, 100*eps)
+%!assert(bessely(alpha,x), yx, 100*eps)
+%!assert(besseli(alpha,x), ix, 100*eps)
+%!assert(besselk(alpha,x), kx, 100*eps)
+%!assert(besselh(alpha,1,x), jx + I*yx, 100*eps)
+%!assert(besselh(alpha,2,x), jx - I*yx, 100*eps)
+%!
+%!assert(besselj(alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
+%!assert(bessely(alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
+%!assert(besseli(alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
+%!assert(besselk(alpha,x,1), kx*exp(x), 100*eps)
+%!assert(besselh(alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
+%!assert(besselh(alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
+%!
+%!assert(besselj(-alpha,x), -jx, 100*eps)
+%!assert(bessely(-alpha,x), -yx, 100*eps)
+%!assert(besseli(-alpha,x), ix, 100*eps)
+%!assert(besselk(-alpha,x), kx, 100*eps)
+%!assert(besselh(-alpha,1,x), -(jx + I*yx), 100*eps)
+%!assert(besselh(-alpha,2,x), -(jx - I*yx), 100*eps)
+%!
+%!assert(besselj(-alpha,x,1), -jx*exp(-abs(imag(x))), 100*eps)
+%!assert(bessely(-alpha,x,1), -yx*exp(-abs(imag(x))), 100*eps)
+%!assert(besseli(-alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
+%!assert(besselk(-alpha,x,1), kx*exp(x), 100*eps)
+%!assert(besselh(-alpha,1,x,1), -(jx + I*yx)*exp(-I*x), 100*eps)
+%!assert(besselh(-alpha,2,x,1), -(jx - I*yx)*exp(I*x), 100*eps)
+%!
+%! x *= -1;
+%! jx = -jx;
+%! yx = 0.7560554967536709968379029772 - 0.4332007820782270495333780070*I;
+%! ix = -ix;
+%! kx = -0.2682271463934492027663765197 - 1.490278591297463775542004240*I;
+%!
+%!assert(besselj(alpha,x), jx, 100*eps)
+%!assert(bessely(alpha,x), yx, 100*eps)
+%!assert(besseli(alpha,x), ix, 100*eps)
+%!assert(besselk(alpha,x), kx, 100*eps)
+%!assert(besselh(alpha,1,x), jx + I*yx, 100*eps)
+%!assert(besselh(alpha,2,x), jx - I*yx, 100*eps)
+%!
+%!assert(besselj(alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
+%!assert(bessely(alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
+%!assert(besseli(alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
+%!assert(besselk(alpha,x,1), kx*exp(x), 100*eps)
+%!assert(besselh(alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
+%!assert(besselh(alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
+%!
+%! # Bessel functions, fractional order, positive and negative x
+%!
+%! alpha = 3.5; x = 2.75;
+%! jx = 0.1691636439842384154644784389;
+%! yx = -0.8301381935499356070267953387;
+%! ix = 0.3930540878794826310979363668;
+%! kx = 0.2844099013460621170288192503;
+%!
+%!assert(besselj(alpha,x), jx, 100*eps)
+%!assert(bessely(alpha,x), yx, 100*eps)
+%!assert(besseli(alpha,x), ix, 100*eps)
+%!assert(besselk(alpha,x), kx, 100*eps)
+%!assert(besselh(alpha,1,x), jx + I*yx, 100*eps)
+%!assert(besselh(alpha,2,x), jx - I*yx, 100*eps)
+%!
+%!assert(besselj(alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
+%!assert(bessely(alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
+%!assert(besseli(alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
+%!assert(besselk(alpha,x,1), kx*exp(x), 100*eps)
+%!assert(besselh(alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
+%!assert(besselh(alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
+%!
+%! nix = 0.2119931212254662995364461998;
+%!
+%!assert(besselj(-alpha,x), yx, 100*eps)
+%!assert(bessely(-alpha,x), -jx, 100*eps)
+%!assert(besseli(-alpha,x), nix, 100*eps)
+%!assert(besselk(-alpha,x), kx, 100*eps)
+%!assert(besselh(-alpha,1,x), -I*(jx + I*yx), 100*eps)
+%!assert(besselh(-alpha,2,x), I*(jx - I*yx), 100*eps)
+%!
+%!assert(besselj(-alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
+%!assert(bessely(-alpha,x,1), -jx*exp(-abs(imag(x))), 100*eps)
+%!assert(besseli(-alpha,x,1), nix*exp(-abs(real(x))), 100*eps)
+%!assert(besselk(-alpha,x,1), kx*exp(x), 100*eps)
+%!assert(besselh(-alpha,1,x,1), -I*(jx + I*yx)*exp(-I*x), 100*eps)
+%!assert(besselh(-alpha,2,x,1), I*(jx - I*yx)*exp(I*x), 100*eps)
+%!
+%! x *= -1;
+%! jx *= -I;
+%! yx = -0.8301381935499356070267953387*I;
+%! ix *= -I;
+%! kx = -0.9504059335995575096509874508*I;
+%!
+%!assert(besselj(alpha,x), jx, 100*eps)
+%!assert(bessely(alpha,x), yx, 100*eps)
+%!assert(besseli(alpha,x), ix, 100*eps)
+%!assert(besselk(alpha,x), kx, 100*eps)
+%!assert(besselh(alpha,1,x), jx + I*yx, 100*eps)
+%!assert(besselh(alpha,2,x), jx - I*yx, 100*eps)
+%!
+%!assert(besselj(alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
+%!assert(bessely(alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
+%!assert(besseli(alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
+%!assert(besselk(alpha,x,1), kx*exp(x), 100*eps)
+%!assert(besselh(alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
+%!assert(besselh(alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
+%!
+%! # Bessel functions, even order, complex x
+%!
+%! alpha = 2; x = 1.25 + 3.625 * I;
+%! jx = -1.299533366810794494030065917 + 4.370833116012278943267479589*I;
+%! yx = -4.370357232383223896393056727 - 1.283083391453582032688834041*I;
+%! ix = -0.6717801680341515541002273932 - 0.2314623443930774099910228553*I;
+%! kx = -0.01108009888623253515463783379 + 0.2245218229358191588208084197*I;
+%!
+%!assert(besselj(alpha,x), jx, 100*eps)
+%!assert(bessely(alpha,x), yx, 100*eps)
+%!assert(besseli(alpha,x), ix, 100*eps)
+%!assert(besselk(alpha,x), kx, 100*eps)
+%!assert(besselh(alpha,1,x), jx + I*yx, 100*eps)
+%!assert(besselh(alpha,2,x), jx - I*yx, 100*eps)
+%!
+%!assert(besselj(alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
+%!assert(bessely(alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
+%!assert(besseli(alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
+%!assert(besselk(alpha,x,1), kx*exp(x), 100*eps)
+%!assert(besselh(alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
+%!assert(besselh(alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
+%!
+%!assert(besselj(-alpha,x), jx, 100*eps)
+%!assert(bessely(-alpha,x), yx, 100*eps)
+%!assert(besseli(-alpha,x), ix, 100*eps)
+%!assert(besselk(-alpha,x), kx, 100*eps)
+%!assert(besselh(-alpha,1,x), jx + I*yx, 100*eps)
+%!assert(besselh(-alpha,2,x), jx - I*yx, 100*eps)
+%!
+%!assert(besselj(-alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
+%!assert(bessely(-alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
+%!assert(besseli(-alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
+%!assert(besselk(-alpha,x,1), kx*exp(x), 100*eps)
+%!assert(besselh(-alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
+%!assert(besselh(-alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
+%!
+%! # Bessel functions, odd order, complex x
+%!
+%! alpha = 3; x = 2.5 + 1.875 * I;
+%! jx = 0.1330721523048277493333458596 + 0.5386295217249660078754395597*I;
+%! yx = -0.6485072392105829901122401551 + 0.2608129289785456797046996987*I;
+%! ix = -0.6182064685486998097516365709 + 0.4677561094683470065767989920*I;
+%! kx = -0.1568585587733540007867882337 - 0.05185853709490846050505141321*I;
+%!
+%!assert(besselj(alpha,x), jx, 100*eps)
+%!assert(bessely(alpha,x), yx, 100*eps)
+%!assert(besseli(alpha,x), ix, 100*eps)
+%!assert(besselk(alpha,x), kx, 100*eps)
+%!assert(besselh(alpha,1,x), jx + I*yx, 100*eps)
+%!assert(besselh(alpha,2,x), jx - I*yx, 100*eps)
+%!
+%!assert(besselj(alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
+%!assert(bessely(alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
+%!assert(besseli(alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
+%!assert(besselk(alpha,x,1), kx*exp(x), 100*eps)
+%!assert(besselh(alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
+%!assert(besselh(alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
+%!
+%!assert(besselj(-alpha,x), -jx, 100*eps)
+%!assert(bessely(-alpha,x), -yx, 100*eps)
+%!assert(besseli(-alpha,x), ix, 100*eps)
+%!assert(besselk(-alpha,x), kx, 100*eps)
+%!assert(besselh(-alpha,1,x), -(jx + I*yx), 100*eps)
+%!assert(besselh(-alpha,2,x), -(jx - I*yx), 100*eps)
+%!
+%!assert(besselj(-alpha,x,1), -jx*exp(-abs(imag(x))), 100*eps)
+%!assert(bessely(-alpha,x,1), -yx*exp(-abs(imag(x))), 100*eps)
+%!assert(besseli(-alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
+%!assert(besselk(-alpha,x,1), kx*exp(x), 100*eps)
+%!assert(besselh(-alpha,1,x,1), -(jx + I*yx)*exp(-I*x), 100*eps)
+%!assert(besselh(-alpha,2,x,1), -(jx - I*yx)*exp(I*x), 100*eps)
+%!
+%! # Bessel functions, fractional order, complex x
+%!
+%! alpha = 3.5; x = 1.75 + 4.125 * I;
+%! jx = -3.018566131370455929707009100 - 0.7585648436793900607704057611*I;
+%! yx = 0.7772278839106298215614791107 - 3.018518722313849782683792010*I;
+%! ix = 0.2100873577220057189038160913 - 0.6551765604618246531254970926*I;
+%! kx = 0.1757147290513239935341488069 + 0.08772348296883849205562558311*I;
+%!
+%!assert(besselj(alpha,x), jx, 100*eps)
+%!assert(bessely(alpha,x), yx, 100*eps)
+%!assert(besseli(alpha,x), ix, 100*eps)
+%!assert(besselk(alpha,x), kx, 100*eps)
+%!assert(besselh(alpha,1,x), jx + I*yx, 100*eps)
+%!assert(besselh(alpha,2,x), jx - I*yx, 100*eps)
+%!
+%!assert(besselj(alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
+%!assert(bessely(alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
+%!assert(besseli(alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
+%!assert(besselk(alpha,x,1), kx*exp(x), 100*eps)
+%!assert(besselh(alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
+%!assert(besselh(alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
+%!
+%! nix = 0.09822388691172060573913739253 - 0.7110230642207380127317227407*I;
+%!
+%!assert(besselj(-alpha,x), yx, 100*eps)
+%!assert(bessely(-alpha,x), -jx, 100*eps)
+%!assert(besseli(-alpha,x), nix, 100*eps)
+%!assert(besselk(-alpha,x), kx, 100*eps)
+%!assert(besselh(-alpha,1,x), -I*(jx + I*yx), 100*eps)
+%!assert(besselh(-alpha,2,x), I*(jx - I*yx), 100*eps)
+%!
+%!assert(besselj(-alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
+%!assert(bessely(-alpha,x,1), -jx*exp(-abs(imag(x))), 100*eps)
+%!assert(besseli(-alpha,x,1), nix*exp(-abs(real(x))), 100*eps)
+%!assert(besselk(-alpha,x,1), kx*exp(x), 100*eps)
+%!assert(besselh(-alpha,1,x,1), -I*(jx + I*yx)*exp(-I*x), 100*eps)
+%!assert(besselh(-alpha,2,x,1), I*(jx - I*yx)*exp(I*x), 100*eps)
+*/
+
+/*
;;; Local Variables: ***
;;; mode: C++ ***
;;; End: ***
-------------- next part --------------
# HG changeset patch
# User Brian Gough <bjg at gnu.org>
# Date 1225198580 14400
# Node ID 6b635febc6d52344064c731e0a2d83715b25af08
# Parent 16699ae4ef0a65eed6127d239ee87239f2dd80df
fix scaling factor for negative alpha in zbesi,cbesi
diff -r 16699ae4ef0a -r 6b635febc6d5 liboctave/ChangeLog
--- a/liboctave/ChangeLog Tue Oct 28 08:54:01 2008 -0400
+++ b/liboctave/ChangeLog Tue Oct 28 08:56:20 2008 -0400
@@ -1,3 +1,8 @@
+2008-10-28 Brian Gough <bjg at gnu.org>
+
+ * lo-specfun.cc (zbesi): Fix scaling factor for negative alpha.
+ (cbesi): Likewise.
+
2008-10-23 John Swensen <jpswensen at comcast.net>
* oct-shlib.cc (octave_dyld_shlib::open): Call NSLinkEditError to
diff -r 16699ae4ef0a -r 6b635febc6d5 liboctave/lo-specfun.cc
--- a/liboctave/lo-specfun.cc Tue Oct 28 08:54:01 2008 -0400
+++ b/liboctave/lo-specfun.cc Tue Oct 28 08:56:20 2008 -0400
@@ -819,8 +819,16 @@
if (ierr == 0 || ierr == 3)
{
- tmp += (2.0 / M_PI) * sin (M_PI * alpha)
+ Complex tmp2 = (2.0 / M_PI) * sin (M_PI * alpha)
* zbesk (z, alpha, kode, ierr);
+
+ if (kode == 2)
+ {
+ /* Compensate for different scaling factor of besk */
+ tmp2 *= exp(-z - std::abs(z.real()));
+ }
+
+ tmp += tmp2;
retval = bessel_return_value (tmp, ierr);
}
@@ -1438,8 +1446,16 @@
if (ierr == 0 || ierr == 3)
{
- tmp += static_cast<float> (2.0 / M_PI) * sinf (static_cast<float> (M_PI) * alpha)
+ FloatComplex tmp2 = static_cast<float> (2.0 / M_PI) * sinf (static_cast<float> (M_PI) * alpha)
* cbesk (z, alpha, kode, ierr);
+
+ if (kode == 2)
+ {
+ /* Compensate for different scaling factor of besk */
+ tmp2 *= exp(-z - std::abs(z.real()));
+ }
+
+ tmp += tmp2;
retval = bessel_return_value (tmp, ierr);
}
More information about the Bug-octave
mailing list