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

scipy-svn@scip... scipy-svn@scip...
Thu Jan 15 18:04:27 CST 2009


Author: ptvirtan
Date: 2009-01-15 18:04:07 -0600 (Thu, 15 Jan 2009)
New Revision: 5465

Modified:
   trunk/scipy/special/specfun/specfun.f
   trunk/scipy/special/tests/test_basic.py
Log:
Fix #679: specfun's STVHV used Y_|v| instead of Y_v => wrong results for v<0

specfun's STVHV uses asymptotic expansion to calculate H_v - Y_v, and
adds Y_v to the result. But it actually added Y_|v| to the result, giving
wrong results for v<0. This commit fixes this by computing Y_v for v<0
via symmetry (integer v) or via a relation between Y_v and Y_|v|, J_|v|.

The Bessel function J_|v|(z) is calculated using forward recurrence,
since this can be done easily. In general this would be numerically
unstable; however, STVHV is restricted to range |v| < 12.5, and since
the asymptotic expansion is triggered only for z > 20, the recurrence
is stable enough. Accuracy 1e-8 appears to be guaranteed.

Modified: trunk/scipy/special/specfun/specfun.f
===================================================================
--- trunk/scipy/special/specfun/specfun.f	2009-01-14 07:29:39 UTC (rev 5464)
+++ trunk/scipy/special/specfun/specfun.f	2009-01-16 00:04:07 UTC (rev 5465)
@@ -12665,6 +12665,7 @@
 C       Input :  v  --- Order of Hv(x)  ( -8.0 ≤ v ≤ 12.5 )
 C                x  --- Argument of Hv(x) ( x ≥ 0 )
 C       Output:  HV --- Hv(x)
+C       Note: numerically unstable away from the above range for `v`
 C       Routine called: GAMMA2 to compute the gamma function
 C       =====================================================
 C
@@ -12685,6 +12686,7 @@
         QU0=0.0D0
         PU0=0.0D0
         IF (X.LE.20.0D0) THEN
+C          Power series for Hv (Abramowitz & Stegun 12.1.3)
            V0=V+1.5D0
            CALL GAMMA2(V0,GA)
            S=2.0D0/(DSQRT(PI)*GA)
@@ -12701,6 +12703,7 @@
 10         CONTINUE
 15         HV=(0.5D0*X)**(V+1.0D0)*S
         ELSE
+C          Asymptotic large |z| expansion for Hv - Yv  (Abm & Stg 12.1.29)
            SA=(0.5D0*X)**(V-1.0)/PI
            V0=V+0.5D0
            CALL GAMMA2(V0,GA)
@@ -12715,6 +12718,8 @@
               S=S+R1*GA/GB
 20         CONTINUE
            S0=SA*S
+
+C          Compute Y_(|v|-N)   (Abm & Stg 9.2.6)
            U=DABS(V)
            N=INT(U)
            U0=U-N
@@ -12745,6 +12750,8 @@
            SR=DSQRT(2.0D0/(PI*X))
            BY0=SR*(PU0*DSIN(T0)+QU0*DCOS(T0))
            BY1=SR*(PU1*DSIN(T1)+QU1*DCOS(T1))
+
+C          Compute Y_|v|   (Abm & Stg 9.1.27)
            BF0=BY0
            BF1=BY1
            DO 40 K=2,N
@@ -12754,6 +12761,37 @@
            IF (N.EQ.0) BYV=BY0
            IF (N.EQ.1) BYV=BY1
            IF (N.GT.1) BYV=BF
+
+C          Compute Y_v  (handle the case v < 0 appropriately)
+           IF (V .LT. 0) THEN
+              IF (U0 .EQ. 0) THEN
+C                Use symmetry (Abm & Stg 9.1.5)
+                 BYV=(-1)**N*BYV
+              ELSE
+C                Use relation between Yv & Jv (Abm & Stg 9.1.6)
+
+C                Compute J_(|v|-N) (Abm & Stg 9.2.5)
+                 BJ0=SR*(PU0*DCOS(T0)-QU0*DSIN(T0))
+                 BJ1=SR*(PU1*DCOS(T1)-QU1*DSIN(T1))
+C                Forward recurrence for J_|v| (Abm & Stg 9.1.27)
+C                It's OK for the limited range -8.0 ≤ v ≤ 12.5, 
+C                since x >= 20 here; but would be unstable for v <~ -20
+                 BF0=BJ0
+                 BF1=BJ1
+                 DO 50 K=2,N
+                    BF=2.0D0*(K-1.0+U0)/X*BF1-BF0
+                    BF0=BF1
+50                  BF1=BF
+                 IF (N.EQ.0) BJV=BJ0
+                 IF (N.EQ.1) BJV=BJ1
+                 IF (N.GT.1) BJV=BF
+
+C                Compute Y_v    (Abm & Stg 9.1.6)
+                 BYV = DCOS(V*PI)*BYV + DSIN(-V*PI)*BJV
+              END IF
+           END IF
+
+C          Compute H_v
            HV=BYV+S0
         ENDIF
         RETURN
@@ -12762,4 +12800,3 @@
 
 
 C       **********************************
-

Modified: trunk/scipy/special/tests/test_basic.py
===================================================================
--- trunk/scipy/special/tests/test_basic.py	2009-01-14 07:29:39 UTC (rev 5464)
+++ trunk/scipy/special/tests/test_basic.py	2009-01-16 00:04:07 UTC (rev 5465)
@@ -24,7 +24,6 @@
 from numpy import array
 
 from numpy.testing import *
-
 from scipy.special import *
 import scipy.special._cephes as cephes
 
@@ -2005,5 +2004,31 @@
         sy3 = sph_yn(1,.2)[1][1]
         assert_almost_equal(sy3,sphpy,4) #compare correct derivative val. (correct =-system val).
 
+class TestStruve(object):
+    def _series(self, v, z, n=100):
+        """Compute Struve function & error estimate from its power series."""
+        k = arange(0, n)
+        r = (-1)**k * (.5*z)**(2*k+v+1)/gamma(k+1.5)/gamma(k+v+1.5)
+        err = abs(r).max() * finfo(float_).eps * n
+        return r.sum(), err
+
+    def test_vs_series(self):
+        """Check Struve function versus its power series"""
+        for v in [-10, -7.99, -3.4, -1, 0, 1, 3.4, 12.49, 16]:
+            for z in [1, 10, 19, 21, 30]:
+                value, err = self._series(v, z)
+                assert allclose(struve(v, z), value, atol=err), (v, z)
+
+    def test_some_values(self):
+        assert_almost_equal(struve(-7.99, 21), 0.0467547614113, decimal=8)
+        assert_almost_equal(struve(-8.01, 21), 0.0398716951023, decimal=9)
+        assert_almost_equal(struve(-3.0, 200), 0.0142134427432, decimal=13)
+
+    def test_regression_679(self):
+        """Regression test for #679"""
+        assert_almost_equal(struve(-1.0, 20 - 1e-8), struve(-1.0, 20 + 1e-8))
+        assert_almost_equal(struve(-2.0, 20 - 1e-8), struve(-2.0, 20 + 1e-8))
+        assert_almost_equal(struve(-4.3, 20 - 1e-8), struve(-4.3, 20 + 1e-8))
+
 if __name__ == "__main__":
     run_module_suite()



More information about the Scipy-svn mailing list