naive quadrature using adaptive subdivision

fredrik-johansson · fredrik-johansson · commit 7123a67a1d40 · 2021-02-26T13:05:54.000+01:00
diff --git a/doc/source/calculus/integration.txt b/doc/source/calculus/integration.txt
@@ -6,6 +6,11 @@ Standard quadrature (``quad``)
 
 .. autofunction:: mpmath.quad
 
+Quadrature with subdivision (``quadsubdiv``)
+............................................
+
+.. autofunction:: mpmath.quadsubdiv
+
 Oscillatory quadrature (``quadosc``)
 ....................................
 
diff --git a/mpmath/__init__.py b/mpmath/__init__.py
@@ -103,6 +103,7 @@
 quadgl = mp.quadgl
 quadts = mp.quadts
 quadosc = mp.quadosc
+quadsubdiv = mp.quadsubdiv
 
 invertlaplace = mp.invertlaplace
 invlaptalbot = mp.invlaptalbot
diff --git a/mpmath/calculus/quadrature.py b/mpmath/calculus/quadrature.py
@@ -718,6 +718,8 @@ def quad(ctx, f, *points, **kwargs):
             >>> quad(f, [-100, 0, 100])   # Also good
             3.12159332021646
 
+        For such integr
+
         **References**
 
         1. http://mathworld.wolfram.com/DoubleIntegral.html
@@ -1006,6 +1008,108 @@ def term(k):
         s += ctx.nsum(term, [n, ctx.inf])
         return s
 
+    def quadsubdiv(ctx, f, interval, tol=None, maxintervals=None, **kwargs):
+        """
+        Computes the integral of *f* over the interval or path specified
+        by *interval*, using :func:`~mpmath.quad` together with adaptive
+        subdivision of the interval.
+
+        This function gives an accurate answer for some integrals where
+        :func:`~mpmath.quad` fails::
+
+            >>> mp.dps = 15; mp.pretty = True
+            >>> quad(lambda x: abs(sin(x)), [0, 2*pi])
+            3.99900894176779
+            >>> quadsubdiv(lambda x: abs(sin(x)), [0, 2*pi])
+            4.0
+            >>> quadsubdiv(sin, [0, 1000])
+            0.437620923709297
+            >>> quadsubdiv(lambda x: 1/(1+x**2), [-100, 100])
+            3.12159332021646
+            >>> quadsubdiv(lambda x: ceil(x), [0, 100])
+            5050.0
+            >>> quadsubdiv(lambda x: sin(x+exp(x)), [0,8])
+            0.347400172657248
+
+        The argument *maxintervals* can be set to limit the permissible
+        subdivision::
+
+            >>> quadsubdiv(lambda x: sin(x**2), [0,100], maxintervals=5, error=True)
+            (-5.40487904307774, 5.011)
+            >>> quadsubdiv(lambda x: sin(x**2), [0,100], maxintervals=100, error=True)
+            (0.631417921866934, 1.10101120134116e-17)
+
+        Subdivision does not guarantee a correct answer since, the error
+        estimate on subintervals may be inaccurate::
+
+            >>> quadsubdiv(lambda x: sech(10*x-2)**2 + sech(100*x-40)**4 + sech(1000*x-600)**6, [0,1], error=True)
+            (0.209736068833883, 1.00011000000001e-18)
+            >>> mp.dps = 20
+            >>> quadsubdiv(lambda x: sech(10*x-2)**2 + sech(100*x-40)**4 + sech(1000*x-600)**6, [0,1], error=True)
+            (0.21080273550054927738, 2.200000001e-24)
+
+        The second answer is correct. We can get an accurate result at lower
+        precision by forcing a finer initial subdivision::
+
+            >>> quadsubdiv(lambda x: sech(10*x-2)**2 + sech(100*x-40)**4 + sech(1000*x-600)**6, linspace(0,1,5))
+            0.210802735500549
+
+        The following integral is too oscillatory for convergence, but we can get a
+        reasonable estimate::
+
+            >>> v, err = fp.quadsubdiv(lambda x: fp.sin(1/x), [0,1], error=True)
+            >>> round(v, 6), round(err, 6)
+            (0.504067, 1e-06)
+            >>> sin(1) - ci(1)
+            0.504067061906928
+
+        """
+        queue = []
+        for i in range(len(interval)-1):
+            queue.append((interval[i], interval[i+1]))
+        total = ctx.zero
+        total_error = ctx.zero
+        if maxintervals is None:
+            maxintervals = 10 * ctx.prec
+        count = 0
+        quad_args = kwargs.copy()
+        quad_args["verbose"] = False
+        quad_args["error"] = True
+        if tol is None:
+            tol = +ctx.eps
+        orig = ctx.prec
+        try:
+            ctx.prec += 5
+            while queue:
+                a, b = queue.pop()
+                s, err = ctx.quad(f, [a, b], **quad_args)
+                if kwargs.get("verbose"):
+                    print("subinterval", count, a, b, err)
+                if err < tol or count > maxintervals:
+                    total += s
+                    total_error += err
+                else:
+                    count += 1
+                    if count == maxintervals and kwargs.get("verbose"):
+                        print("warning: number of intervals exceeded maxintervals")
+                    if a == -ctx.inf and b == ctx.inf:
+                        m = 0
+                    elif a == -ctx.inf:
+                        m = min(b-1, 2*b)
+                    elif b == ctx.inf:
+                        m = max(a+1, 2*a)
+                    else:
+                        m = a + (b - a) / 2
+                    queue.append((a, m))
+                    queue.append((m, b))
+        finally:
+            ctx.prec = orig
+        if kwargs.get("error"):
+            return +total, +total_error
+        else:
+            return +total
+
 if __name__ == '__main__':
     import doctest
     doctest.testmod()
+
diff --git a/mpmath/functions/elliptic.py b/mpmath/functions/elliptic.py
@@ -586,8 +586,8 @@ def RJ_calc(ctx, x, y, z, p, r, integration):
             N += margin
             F = lambda t: 1/(ctx.sqrt(t+x)*ctx.sqrt(t+y)*ctx.sqrt(t+z)*(t+p))
             if integration == 2:
-                return 1.5 * ctx.quad(F, [0, N, ctx.inf])
-            initial_integral = 1.5 * ctx.quad(F, [0, N])
+                return 1.5 * ctx.quadsubdiv(F, [0, N, ctx.inf])
+            initial_integral = 1.5 * ctx.quadsubdiv(F, [0, N])
             x += N; y += N; z += N; p += N
     xm,ym,zm,pm = x,y,z,p
     A0 = Am = (x + y + z + 2*p)/5
diff --git a/mpmath/tests/test_elliptic.py b/mpmath/tests/test_elliptic.py
@@ -654,9 +654,11 @@ def test_elliptic_integrals():
     assert elliprj(0.3068, -4.037+0.0632j, 1.654, -0.9609).ae(-1.17157627949475577 - 0.069182614173988811j, abs_eps=1e-8)
     assert elliprj(0.3068, -4.037+0.00632j, 1.654, -0.9609).ae(-1.17337595670549633 - 0.0623069224526925j, abs_eps=1e-8)
 
-    # these don't work because of inaccurate integration
-    # assert elliprj(0.3068, -4.037-0.0632j, 1.654, -0.9609).ae(1.77940452391261626 + 0.0388711305592447234j, abs_eps=1e-8)
-    # assert elliprj(0.3068, -4.037-0.00632j, 1.654, -0.9609).ae(1.77806722756403055 + 0.0592749824572262329j, abs_eps=1e-8)
+    # these require accurate integration
+    assert elliprj(0.3068, -4.037-0.0632j, 1.654, -0.9609).ae(1.77940452391261626 + 0.0388711305592447234j)
+    assert elliprj(0.3068, -4.037-0.00632j, 1.654, -0.9609).ae(1.77806722756403055 + 0.0592749824572262329j)
+    # issue #571
+    assert ellippi(2.1 + 0.94j, 2.3 + 0.98j, 2.5 + 0.01j).ae(-0.40652414240811963438 + 2.1547659461404749309j)
 
     assert ellippi(2.0-1.0j, 2.0+1.0j).ae(1.8578723151271115 - 1.18642180609983531j)
     assert ellippi(2.0-0.5j, 0.5+1.0j).ae(0.936761970766645807 - 1.61876787838890786j)
