[Scipy-svn] r5508 - in trunk/scipy/special: . specfun tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Mon Jan 19 13:49:36 EST 2009
Author: ptvirtan
Date: 2009-01-19 12:48:55 -0600 (Mon, 19 Jan 2009)
New Revision: 5508
Modified:
trunk/scipy/special/basic.py
trunk/scipy/special/specfun.pyf
trunk/scipy/special/specfun/specfun.f
trunk/scipy/special/tests/test_basic.py
Log:
Fixing #852: make specfun.f:JYZO work also for N >= 101.
In addition to removing the checks for N, this requires rewriting
the Newton iterations in JYZO to be more robust, because for large
N the initial guess used between successive zeros becomes inaccurate.
Additional tests for jn/yn zeros are included.
Modified: trunk/scipy/special/basic.py
===================================================================
--- trunk/scipy/special/basic.py 2009-01-19 18:48:21 UTC (rev 5507)
+++ trunk/scipy/special/basic.py 2009-01-19 18:48:55 UTC (rev 5508)
@@ -81,7 +81,7 @@
raise ValueError, "Arguments must be integers."
if (nt <=0):
raise ValueError, "nt > 0"
- return specfun.jyzo(n,nt)
+ return specfun.jyzo(abs(n),nt)
def jn_zeros(n,nt):
"""Compute nt zeros of the Bessel function Jn(x).
Modified: trunk/scipy/special/specfun/specfun.f
===================================================================
--- trunk/scipy/special/specfun/specfun.f 2009-01-19 18:48:21 UTC (rev 5507)
+++ trunk/scipy/special/specfun/specfun.f 2009-01-19 18:48:55 UTC (rev 5508)
@@ -10408,7 +10408,7 @@
C ======================================================
C Purpose: Compute the zeros of Bessel functions Jn(x),
C Yn(x), and their derivatives
-C Input : n --- Order of Bessel functions ( n ≤ 101 )
+C Input : n --- Order of Bessel functions (n >= 0)
C NT --- Number of zeros (roots)
C Output: RJ0(L) --- L-th zero of Jn(x), L=1,2,...,NT
C RJ1(L) --- L-th zero of Jn'(x), L=1,2,...,NT
@@ -10431,14 +10431,23 @@
ENDIF
L=0
C 2) iterate
+ XGUESS=X
10 X0=X
CALL JYNDD(N,X,BJN,DJN,FJN,BYN,DYN,FYN)
X=X-BJN/DJN
+ IF (X-X0.LT.-1) X=X0-1
+ IF (X-X0.GT.1) X=X0+1
IF (DABS(X-X0).GT.1.0D-11) GO TO 10
C 3) initial guess for j_{N,L+1}
+ IF (L.GE.1 .AND. X.LE.RJ0(L)+0.5) THEN
+ X=XGUESS+PI
+ XGUESS=X
+ GO TO 10
+ END IF
L=L+1
RJ0(L)=X
- X=X+PI+(0.0972+0.0679*N-0.000354*N**2)/L
+C XXX: should have a better initial guess for large N ~> 100 here
+ X=X+PI+MAX((0.0972+0.0679*N-0.000354*N**2)/L, 0d0)
IF (L.LT.NT) GO TO 10
C -- Newton method for j_{N,L}'
IF (N.LE.20) THEN
@@ -10448,43 +10457,68 @@
ENDIF
IF (N.EQ.0) X=3.8317
L=0
+ XGUESS=X
15 X0=X
CALL JYNDD(N,X,BJN,DJN,FJN,BYN,DYN,FYN)
X=X-DJN/FJN
+ IF (X-X0.LT.-1) X=X0-1
+ IF (X-X0.GT.1) X=X0+1
IF (DABS(X-X0).GT.1.0D-11) GO TO 15
+ IF (L.GE.1 .AND. X.LE.RJ1(L)+0.5) THEN
+ X=XGUESS+PI
+ XGUESS=X
+ GO TO 15
+ END IF
L=L+1
RJ1(L)=X
- X=X+PI+(0.4955+0.0915*N-0.000435*N**2)/L
+C XXX: should have a better initial guess for large N ~> 100 here
+ X=X+PI+MAX((0.4955+0.0915*N-0.000435*N**2)/L,0d0)
IF (L.LT.NT) GO TO 15
+C -- Newton method for y_{N,L}
IF (N.LE.20) THEN
X=1.19477+1.08933*N
ELSE
X=N+0.93158*N**0.33333+0.26035/N**0.33333
ENDIF
-C -- Newton method for y_{N,L}
L=0
+ XGUESS=X
20 X0=X
CALL JYNDD(N,X,BJN,DJN,FJN,BYN,DYN,FYN)
X=X-BYN/DYN
+ IF (X-X0.LT.-1) X=X0-1
+ IF (X-X0.GT.1) X=X0+1
IF (DABS(X-X0).GT.1.0D-11) GO TO 20
+ IF (L.GE.1 .AND. X.LE.RY0(L)+0.5) THEN
+ X=XGUESS+PI
+ XGUESS=X
+ GO TO 20
+ END IF
L=L+1
RY0(L)=X
- X=X+PI+(0.312+0.0852*N-0.000403*N**2)/L
+C XXX: should have a better initial guess for large N ~> 100 here
+ X=X+PI+MAX((0.312+0.0852*N-0.000403*N**2)/L,0d0)
IF (L.LT.NT) GO TO 20
+C -- Newton method for y_{N,L}'
IF (N.LE.20) THEN
X=2.67257+1.16099*N
ELSE
X=N+1.8211*N**0.33333+0.94001/N**0.33333
ENDIF
-C -- Newton method for y_{N,L}'
L=0
+ XGUESS=X
25 X0=X
CALL JYNDD(N,X,BJN,DJN,FJN,BYN,DYN,FYN)
X=X-DYN/FYN
IF (DABS(X-X0).GT.1.0D-11) GO TO 25
+ IF (L.GE.1 .AND. X.LE.RY1(L)+0.5) THEN
+ X=XGUESS+PI
+ XGUESS=X
+ GO TO 25
+ END IF
L=L+1
RY1(L)=X
- X=X+PI+(0.197+0.0643*N-0.000286*N**2)/L
+C XXX: should have a better initial guess for large N ~> 100 here
+ X=X+PI+MAX((0.197+0.0643*N-0.000286*N**2)/L,0d0)
IF (L.LT.NT) GO TO 25
RETURN
END
Modified: trunk/scipy/special/specfun.pyf
===================================================================
--- trunk/scipy/special/specfun.pyf 2009-01-19 18:48:21 UTC (rev 5507)
+++ trunk/scipy/special/specfun.pyf 2009-01-19 18:48:55 UTC (rev 5508)
@@ -410,7 +410,7 @@
! rswfo
! ch12n
subroutine jyzo(n,nt,rj0,rj1,ry0,ry1) ! in :specfun:specfun.f
- integer intent(in), check(abs(n)<=101) :: n
+ integer intent(in), check(n>=0) :: n
integer intent(in), check(nt>0) :: nt
double precision intent(out),dimension(nt),depend(nt) :: rj0
double precision intent(out),dimension(nt),depend(nt) :: rj1
Modified: trunk/scipy/special/tests/test_basic.py
===================================================================
--- trunk/scipy/special/tests/test_basic.py 2009-01-19 18:48:21 UTC (rev 5507)
+++ trunk/scipy/special/tests/test_basic.py 2009-01-19 18:48:55 UTC (rev 5508)
@@ -1399,16 +1399,38 @@
13.32369,
16.47063]),4)
- def test_jn_zeros_large(self):
+ jn102 = jn_zeros(102,5)
+ assert_tol_equal(jn102, array([110.89174935992040343,
+ 117.83464175788308398,
+ 123.70194191713507279,
+ 129.02417238949092824,
+ 134.00114761868422559]), rtol=1e-13)
+
+ jn301 = jn_zeros(301,5)
+ assert_tol_equal(jn301, array([313.59097866698830153,
+ 323.21549776096288280,
+ 331.22338738656748796,
+ 338.39676338872084500,
+ 345.03284233056064157]), rtol=1e-13)
+
+ def test_jn_zeros_slow(self):
jn0 = jn_zeros(0, 300)
assert_tol_equal(jn0[260-1], 816.02884495068867280, rtol=1e-13)
assert_tol_equal(jn0[280-1], 878.86068707124422606, rtol=1e-13)
assert_tol_equal(jn0[300-1], 941.69253065317954064, rtol=1e-13)
+
jn10 = jn_zeros(10, 300)
assert_tol_equal(jn10[260-1], 831.67668514305631151, rtol=1e-13)
assert_tol_equal(jn10[280-1], 894.51275095371316931, rtol=1e-13)
assert_tol_equal(jn10[300-1], 957.34826370866539775, rtol=1e-13)
+ jn3010 = jn_zeros(3010,5)
+ assert_tol_equal(jn3010, array([3036.86590780927,
+ 3057.06598526482,
+ 3073.66360690272,
+ 3088.37736494778,
+ 3101.86438139042]), rtol=1e-8)
+
def test_jnjnp_zeros(self):
pass
#jnjp = jnjnp(3)
@@ -1423,6 +1445,8 @@
8.53632,
11.70600,
14.86359]),4)
+ jnp = jnp_zeros(443,5)
+ assert_tol_equal(jvp(443, jnp), 0, atol=1e-15)
def test_jnyn_zeros(self):
jnz = jnyn_zeros(1,5)
@@ -1445,7 +1469,7 @@
6.94150,
10.12340,
13.28576,
- 16.44006])),4)
+ 16.44006])),5)
def test_jvp(self):
jvprim = jvp(2,2)
@@ -1568,10 +1592,16 @@
def test_yn_zeros(self):
an = yn_zeros(4,2)
assert_array_almost_equal(an,array([ 5.64515, 9.36162]),5)
+ an = yn_zeros(443,5)
+ assert_tol_equal(an, [450.13573091578090314, 463.05692376675001542,
+ 472.80651546418663566, 481.27353184725625838,
+ 488.98055964441374646], rtol=1e-19)
def test_ynp_zeros(self):
ao = ynp_zeros(0,2)
assert_array_almost_equal(ao,array([ 2.19714133, 5.42968104]),6)
+ ao = ynp_zeros(443,5)
+ assert_tol_equal(yvp(443, ao), 0, atol=1e-15)
def test_yn(self):
yn2n = yn(1,.2)
More information about the Scipy-svn
mailing list