[Scipy-svn] r5508 - in trunk/scipy/special: . specfun tests

scipy-svn@scip... scipy-svn@scip...
Mon Jan 19 12:49:36 CST 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