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

scipy-svn at scipy.org scipy-svn at scipy.org
Mon Jan 19 13:48:39 EST 2009


Author: ptvirtan
Date: 2009-01-19 12:48:21 -0600 (Mon, 19 Jan 2009)
New Revision: 5507

Modified:
   trunk/scipy/special/specfun/specfun.f
   trunk/scipy/special/tests/test_basic.py
Log:
Fix #852: improve accuracy of jn_zeros by refactoring Specfun JYNDD, JYNB

Code to calculate the value of Bessel functions of integer order was
duplicated in specfun.f routines JYNDD and JYNB. The former routine
didn't do asymptotic expansion, and its results were bad for large
values of x, resulting to errors in `scipy.special.jn_zeros`.

This commit refactors the common code of these two functions to
JYNBH, and makes JYZO Newton tolerances stricter. The changes to JYNDD
also remove the limitation n <= 101.

Improved tests for jn_zeros are included.

Modified: trunk/scipy/special/specfun/specfun.f
===================================================================
--- trunk/scipy/special/specfun/specfun.f	2009-01-19 10:29:13 UTC (rev 5506)
+++ trunk/scipy/special/specfun/specfun.f	2009-01-19 18:48:21 UTC (rev 5507)
@@ -1863,55 +1863,24 @@
 C                BYN ---  Yn(x)
 C                DYN ---  Yn'(x)
 C                FYN ---  Yn"(x)
+C       Routines called:
+C                JYNBH to compute Jn and Yn
 C       ===========================================================
 C
         IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-        DIMENSION BJ(102),BY(102)
-        DO 10 NT=1,900
-           MT=INT(0.5*LOG10(6.28*NT)-NT*LOG10(1.36*DABS(X)/NT))
-           IF (MT.GT.20) GO TO 15
-10      CONTINUE
-15      M=NT
-        BS=0.0D0
-        F0=0.0D0
-        F1=1.0D-35
-        SU=0.0D0
-        F=0.0D0
-        DO 20 K=M,0,-1
-           F=2.0D0*(K+1.0D0)*F1/X-F0
-           IF (K.LE.N+1) BJ(K+1)=F
-           IF (K.EQ.2*INT(K/2)) THEN
-              BS=BS+2.0D0*F
-              IF (K.NE.0) SU=SU+(-1)**(K/2)*F/K
-           ENDIF
-           F0=F1
-20         F1=F
-        DO 25 K=0,N+1
-25         BJ(K+1)=BJ(K+1)/(BS-F)
-        BJN=BJ(N+1)
-        EC=0.5772156649015329D0
-        E0=0.3183098861837907D0
-        S1=2.0D0*E0*(DLOG(X/2.0D0)+EC)*BJ(1)
-        F0=S1-8.0D0*E0*SU/(BS-F)
-        F1=(BJ(2)*F0-2.0D0*E0/X)/BJ(1)
-        BY(1)=F0
-        BY(2)=F1
-        DO 30 K=2,N+1
-           F=2.0D0*(K-1.0D0)*F1/X-F0
-           BY(K+1)=F
-           F0=F1
-30         F1=F
-        BYN=BY(N+1)
-        DJN=-BJ(N+2)+N*BJ(N+1)/X
-        DYN=-BY(N+2)+N*BY(N+1)/X
+        DIMENSION BJ(2),BY(2)
+        CALL JYNBH(N+1,N,X,NM,BJ,BY)
+C       Compute derivatives by differentiation formulas
+        BJN=BJ(1)
+        BYN=BY(1)
+        DJN=-BJ(2)+N*BJ(1)/X
+        DYN=-BY(2)+N*BY(1)/X
         FJN=(N*N/(X*X)-1.0D0)*BJN-DJN/X
         FYN=(N*N/(X*X)-1.0D0)*BYN-DYN/X
         RETURN
         END
 
 
- 
-
 C       **********************************
 
         SUBROUTINE GAM0 (X,GA)
@@ -4253,27 +4222,61 @@
 C                DY(n) --- Yn'(x)
 C                NM --- Highest order computed
 C       Routines called:
+C                JYNBH to calculate the Jn and Yn
+C       =====================================================
+C
+        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+        DIMENSION BJ(0:N),DJ(0:N),BY(0:N),DY(0:N)
+        CALL JYNBH(N,0,X,NM,BJ,BY)
+C       Compute derivatives by differentiation formulas
+        IF (X.LT.1.0D-100) THEN
+           DO 10 K=0,N
+              DJ(K) = 0.0D0
+ 10           DY(K) = 1.0D+300
+           DJ(1)=0.5D0
+        ELSE   
+           DJ(0)=-BJ(1)
+           DO 40 K=1,NM
+ 40           DJ(K)=BJ(K-1)-K/X*BJ(K)
+           DY(0)=-BY(1)
+           DO 50 K=1,NM
+ 50           DY(K)=BY(K-1)-K*BY(K)/X
+        END IF
+        RETURN
+        END
+
+
+C       **********************************
+
+        SUBROUTINE JYNBH(N,NMIN,X,NM,BJ,BY)
+C
+C       =====================================================
+C       Purpose: Compute Bessel functions Jn(x), Yn(x)
+C       Input :  x --- Argument of Jn(x) and Yn(x) ( x ≥ 0 )
+C                n --- Highest order of Jn(x) and Yn(x) computed  ( n ≥ 0 )
+C                nmin -- Lowest order computed  ( nmin ≥ 0 )
+C       Output:  BJ(n-NMIN) --- Jn(x)   ; if indexing starts at 0
+C                BY(n-NMIN) --- Yn(x)   ; if indexing starts at 0
+C                NM --- Highest order computed
+C       Routines called:
 C                MSTA1 and MSTA2 to calculate the starting 
 C                point for backward recurrence
 C       =====================================================
 C
         IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-        DIMENSION BJ(0:N),DJ(0:N),BY(0:N),DY(0:N),
-     &            A(4),B(4),A1(4),B1(4)
+        DIMENSION BJ(0:N-NMIN),BY(0:N-NMIN),A(4),B(4),A1(4),B1(4)
         PI=3.141592653589793D0
         R2P=.63661977236758D0
         NM=N
         IF (X.LT.1.0D-100) THEN
-           DO 10 K=0,N
-              BJ(K)=0.0D0
-              DJ(K)=0.0D0
-              BY(K)=-1.0D+300
-10            DY(K)=1.0D+300
-           BJ(0)=1.0D0
-           DJ(1)=0.5D0
+           DO 10 K=NMIN,N
+              BJ(K-NMIN)=0.0D0
+10            BY(K-NMIN)=-1.0D+300
+           IF (NMIN.EQ.0) BJ(0)=1.0D0
            RETURN
         ENDIF
         IF (X.LE.300.0.OR.N.GT.INT(0.9*X)) THEN
+C          Backward recurrence for Jn
            IF (N.EQ.0) NM=1
            M=MSTA1(X,200)
            IF (M.LT.NM) THEN
@@ -4289,7 +4292,7 @@
            F=0.0D0
            DO 15 K=M,0,-1
               F=2.0D0*(K+1.0D0)/X*F1-F2
-              IF (K.LE.NM) BJ(K)=F
+              IF (K.LE.NM .AND. K.GE.NMIN) BJ(K-NMIN)=F
               IF (K.EQ.2*INT(K/2).AND.K.NE.0) THEN
                  BS=BS+2.0D0*F
                  SU=SU+(-1)**(K/2)*F/K
@@ -4299,14 +4302,19 @@
               F2=F1
 15            F1=F
            S0=BS+F
-           DO 20 K=0,NM
-20            BJ(K)=BJ(K)/S0
+           DO 20 K=NMIN,NM
+20            BJ(K-NMIN)=BJ(K-NMIN)/S0
+C          Estimates for Yn at start of recurrence
+           BJ0 = F1 / S0
+           BJ1 = F2 / S0
            EC=DLOG(X/2.0D0)+0.5772156649015329D0
-           BY0=R2P*(EC*BJ(0)-4.0D0*SU/S0)
-           BY(0)=BY0
-           BY1=R2P*((EC-1.0D0)*BJ(1)-BJ(0)/X-4.0D0*SV/S0)
-           BY(1)=BY1
+           BY0=R2P*(EC*BJ0-4.0D0*SU/S0)
+           BY1=R2P*((EC-1.0D0)*BJ1-BJ0/X-4.0D0*SV/S0)
+           IF (0.GE.NMIN) BY(0-NMIN)=BY0
+           IF (1.GE.NMIN) BY(1-NMIN)=BY1
+           KY=2
         ELSE
+C          Hankel expansion
            DATA A/-.7031250000000000D-01,.1121520996093750D+00,
      &            -.5725014209747314D+00,.6074042001273483D+01/
            DATA B/ .7324218750000000D-01,-.2271080017089844D+00,
@@ -4324,8 +4332,8 @@
            CU=DSQRT(R2P/X)
            BJ0=CU*(P0*DCOS(T1)-Q0*DSIN(T1))
            BY0=CU*(P0*DSIN(T1)+Q0*DCOS(T1))
-           BJ(0)=BJ0
-           BY(0)=BY0
+           IF (0.GE.NMIN) BJ(0-NMIN)=BJ0
+           IF (0.GE.NMIN) BY(0-NMIN)=BY0
            T2=X-0.75D0*PI
            P1=1.0D0
            Q1=0.375D0/X
@@ -4334,25 +4342,21 @@
 30            Q1=Q1+B1(K)*X**(-2*K-1)
            BJ1=CU*(P1*DCOS(T2)-Q1*DSIN(T2))
            BY1=CU*(P1*DSIN(T2)+Q1*DCOS(T2))
-           BJ(1)=BJ1
-           BY(1)=BY1
+           IF (1.GE.NMIN) BJ(1-NMIN)=BJ1
+           IF (1.GE.NMIN) BY(1-NMIN)=BY1
            DO 35 K=2,NM
               BJK=2.0D0*(K-1.0D0)/X*BJ1-BJ0
-              BJ(K)=BJK
+              IF (K.GE.NMIN) BJ(K-NMIN)=BJK
               BJ0=BJ1
 35            BJ1=BJK
+           KY=2
         ENDIF
-        DJ(0)=-BJ(1)
-        DO 40 K=1,NM
-40         DJ(K)=BJ(K-1)-K/X*BJ(K)
-        DO 45 K=2,NM
+C       Forward recurrence for Yn
+        DO 45 K=KY,NM
            BYK=2.0D0*(K-1.0D0)*BY1/X-BY0
-           BY(K)=BYK
+           IF (K.GE.NMIN) BY(K-NMIN)=BYK
            BY0=BY1
 45         BY1=BYK
-        DY(0)=-BY(1)
-        DO 50 K=1,NM
-50         DY(K)=BY(K-1)-K*BY(K)/X
         RETURN
         END
 
@@ -10416,20 +10420,27 @@
 C
         IMPLICIT DOUBLE PRECISION (A-H,O-Z)
         DIMENSION RJ0(NT),RJ1(NT),RY0(NT),RY1(NT)
+        PI=3.141592653589793D0
+C       -- Newton method for j_{N,L}
+C       1) initial guess for j_{N,1}
         IF (N.LE.20) THEN
            X=2.82141+1.15859*N
         ELSE
+C          Abr & Stg (9.5.14)
            X=N+1.85576*N**0.33333+1.03315/N**0.33333
         ENDIF
         L=0
+C       2) iterate
 10      X0=X
         CALL JYNDD(N,X,BJN,DJN,FJN,BYN,DYN,FYN)
         X=X-BJN/DJN
-        IF (DABS(X-X0).GT.1.0D-9) GO TO 10
+        IF (DABS(X-X0).GT.1.0D-11) GO TO 10
+C       3) initial guess for j_{N,L+1}
         L=L+1
         RJ0(L)=X
-        X=X+3.1416+(0.0972+0.0679*N-0.000354*N**2)/L
+        X=X+PI+(0.0972+0.0679*N-0.000354*N**2)/L
         IF (L.LT.NT) GO TO 10
+C       -- Newton method for j_{N,L}'
         IF (N.LE.20) THEN
            X=0.961587+1.07703*N
         ELSE
@@ -10440,38 +10451,40 @@
 15      X0=X
         CALL JYNDD(N,X,BJN,DJN,FJN,BYN,DYN,FYN)
         X=X-DJN/FJN
-        IF (DABS(X-X0).GT.1.0D-9) GO TO 15
+        IF (DABS(X-X0).GT.1.0D-11) GO TO 15
         L=L+1
         RJ1(L)=X
-        X=X+3.1416+(0.4955+0.0915*N-0.000435*N**2)/L
+        X=X+PI+(0.4955+0.0915*N-0.000435*N**2)/L
         IF (L.LT.NT) GO TO 15
         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
 20      X0=X
         CALL JYNDD(N,X,BJN,DJN,FJN,BYN,DYN,FYN)
         X=X-BYN/DYN
-        IF (DABS(X-X0).GT.1.0D-9) GO TO 20
+        IF (DABS(X-X0).GT.1.0D-11) GO TO 20
         L=L+1
         RY0(L)=X
-        X=X+3.1416+(0.312+0.0852*N-0.000403*N**2)/L
+        X=X+PI+(0.312+0.0852*N-0.000403*N**2)/L
         IF (L.LT.NT) GO TO 20
         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
 25      X0=X
         CALL JYNDD(N,X,BJN,DJN,FJN,BYN,DYN,FYN)
         X=X-DYN/FYN
-        IF (DABS(X-X0).GT.1.0D-9) GO TO 25
+        IF (DABS(X-X0).GT.1.0D-11) GO TO 25
         L=L+1
         RY1(L)=X
-        X=X+3.1416+(0.197+0.0643*N-0.000286*N**2)/L 
+        X=X+PI+(0.197+0.0643*N-0.000286*N**2)/L 
         IF (L.LT.NT) GO TO 25
         RETURN
         END

Modified: trunk/scipy/special/tests/test_basic.py
===================================================================
--- trunk/scipy/special/tests/test_basic.py	2009-01-19 10:29:13 UTC (rev 5506)
+++ trunk/scipy/special/tests/test_basic.py	2009-01-19 18:48:21 UTC (rev 5507)
@@ -1399,6 +1399,16 @@
                                               13.32369,
                                               16.47063]),4)
 
+    def test_jn_zeros_large(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)
+
     def test_jnjnp_zeros(self):
         pass
         #jnjp = jnjnp(3)




More information about the Scipy-svn mailing list