Repository: sympy/sympy
Pull request: Symbolic Riemann sums for definite integrals
#### References to other Issues or PRs

Fixes #12707


#### Brief description of what is fixed or changed

Adds a Boolean flag `evaluate` to `as_sum` method, True by default for backward compatibility.

Also allows a symbolic number of subintervals to be used, which is both natural and more interesting for a symbolic library: for some simple functions the sum can be evaluated with a general n, showing that it converges to the integral (and at different rates for different rules).

Concerning the values at endpoints: by the definition of a Riemann sum, direct evaluation `f(a)` or `f(b)` is performed, not a limit as x->a or x->b. For example,
```
Integral(sign(x), (x, 0, 1)).as_sum(1, 'trapezoid')
```
should be (f(0)+f(1))/2 = 1/2, not 1 as it was returned previously.

#### Other comments

A link to Wikipedia is updated, as the previously referenced article was merged.

Previously, passing oo as the number of subintervals or "qwerty" as the method would raise NotImplementedError. Now it raises ValueError, as is appropriate for absurd/invalid inputs. 

Problem statement:
Integral.as_sum() should output a Sum() object
Currently, Integral.as_sum() outputs an evaluated summation instead of an unevaluated expression:

```
In [1]: import sympy as sm

In [2]: t, t0, tf = sm.symbols('t, t0, tf')

In [3]: x = sm.Function('x')(t)

In [4]: y = sm.Function('y')(t)

In [5]: J = sm.Integral((x - y)**2, (t, t0, tf))

In [6]: J.as_sum(20, 'trapezoid')
Out[6]: (-t0/20 + tf/20)*((x(t0/20 + 19*tf/20) - y(t0/20 + 19*tf/20))**2 + (x(t0/10 + 9*tf/10) - y(t0/10 + 9*tf/10))**2 + (x(3*t0/20 + 17*tf/20) - y(3*t0/20 + 17*tf/20))**2 + (x(t0/5 + 4*tf/5) - y(t0/5 + 4*tf/5))**2 + (x(t0/4 + 3*tf/4) - y(t0/4 + 3*tf/4))**2 + (x(3*t0/10 + 7*tf/10) - y(3*t0/10 + 7*tf/10))**2 + (x(7*t0/20 + 13*tf/20) - y(7*t0/20 + 13*tf/20))**2 + (x(2*t0/5 + 3*tf/5) - y(2*t0/5 + 3*tf/5))**2 + (x(9*t0/20 + 11*tf/20) - y(9*t0/20 + 11*tf/20))**2 + (x(t0/2 + tf/2) - y(t0/2 + tf/2))**2 + (x(11*t0/20 + 9*tf/20) - y(11*t0/20 + 9*tf/20))**2 + (x(3*t0/5 + 2*tf/5) - y(3*t0/5 + 2*tf/5))**2 + (x(13*t0/20 + 7*tf/20) - y(13*t0/20 + 7*tf/20))**2 + (x(7*t0/10 + 3*tf/10) - y(7*t0/10 + 3*tf/10))**2 + (x(3*t0/4 + tf/4) - y(3*t0/4 + tf/4))**2 + (x(4*t0/5 + tf/5) - y(4*t0/5 + tf/5))**2 + (x(17*t0/20 + 3*tf/20) - y(17*t0/20 + 3*tf/20))**2 + (x(9*t0/10 + tf/10) - y(9*t0/10 + tf/10))**2 + (x(19*t0/20 + tf/20) - y(19*t0/20 + tf/20))**2 + x(t0)**2/2 - x(t0)*y(t0) + x(tf)**2/2 - x(tf)*y(tf) + y(t0)**2/2 + y(tf)**2/2)
```

For large n this takes a long time to compute. It seems like this should output an unevaluated sum and if the user wants to expand the sum they'd call `.doit()` on the result. It may not be worth deprecating this behavior, but maybe we need to have a `as_unevaluated_sum()` method.

Hints from issue discussion:
Hi @moorepants ! I want to work on this issue .
Go for it!
@moorepants Any kind of help will be appreciated.I'm new to this organization.
Start here: https://github.com/sympy/sympy/wiki/introduction-to-contributing

Read the materials, setup your dev environment, and do the tutorial.
@moorepants I've read the materials.
Coming back to the issue tell me about the changes need to be done and point me to the files to be worked with.
Here is the associated code: https://github.com/sympy/sympy/blob/master/sympy/integrals/integrals.py#L993

I suggest seeing if the reviewers would accept a deprecation so that we can change the behavior of `as_sum` but if not you'll need to create a new method. It would also be worth determining if the `Sum` object has duplicate code for the expansions in it's `doit()` method.

Commit message: Symbolic Riemann sums for definite integrals

Fixes #12707

Adds a Boolean flag `evaluate` to `as_sum` method, True by default for
backward compatibility.

Also allows a symbolic number of subintervals to be used, which is both
natural and more interesting for a symbolic library: for some simple
functions the sum can be evaluated with a general n, showing that it
converges to the integral (and at different rates for different rules).

Concerning the values at endpoints: by the definition of a Riemann sum,
direct evaluation `f(a)` or `f(b)` is performed, not a limit as x->a or
x->b. For example,
```
Integral(sign(x), (x, 0, 1)).as_sum(1, 'trapezoid')
```
should be (f(0)+f(1))/2 = 1/2, not 1 as it was returned previously.

A link to Wikipedia is updated, as the previously referenced article was
merged.

Patch to review:
diff --git a/sympy/integrals/integrals.py b/sympy/integrals/integrals.py
index ea787ea94b80..3e8742828aa6 100644
--- a/sympy/integrals/integrals.py
+++ b/sympy/integrals/integrals.py
@@ -1047,23 +1047,29 @@ def _eval_as_leading_term(self, x):
                 break
         return integrate(leading_term, *self.args[1:])
 
-    def as_sum(self, n, method="midpoint"):
+    def as_sum(self, n, method="midpoint", evaluate=True):
         """
-        Approximates the definite integral by a sum.
+        Approximates a definite integral by a sum.
 
-        method ... one of: left, right, midpoint, trapezoid
+        Arguments
+        ---------
+        n
+            The number of subintervals to use.
+        method
+            One of: 'left', 'right', 'midpoint', 'trapezoid'
+        evaluate
+            If False, returns an unevaluated Sum expression. The default
+            is True, evaluate the sum.
 
-        These are all basically the rectangle method [1], the only difference
-        is where the function value is taken in each interval to define the
-        rectangle.
+        These methods of approximate integration are described in [1].
 
-        [1] http://en.wikipedia.org/wiki/Rectangle_method
+        [1] https://en.wikipedia.org/wiki/Riemann_sum#Methods
 
         Examples
         ========
 
         >>> from sympy import sin, sqrt
-        >>> from sympy.abc import x
+        >>> from sympy.abc import x, n
         >>> from sympy.integrals import Integral
         >>> e = Integral(sin(x), (x, 3, 7))
         >>> e
@@ -1098,9 +1104,8 @@ def as_sum(self, n, method="midpoint"):
         >>> (e.as_sum(2, 'left') + e.as_sum(2, 'right'))/2 == _
         True
 
-        All but the trapexoid method may be used when dealing with a function
-        with a discontinuity. Here, the discontinuity at x = 0 can be avoided
-        by using the midpoint or right-hand method:
+        Here, the discontinuity at x = 0 can be avoided by using the
+        midpoint or right-hand method:
 
         >>> e = Integral(1/sqrt(x), (x, 0, 1))
         >>> e.as_sum(5).n(4)
@@ -1111,12 +1116,24 @@ def as_sum(self, n, method="midpoint"):
         2.000
 
         The left- or trapezoid method will encounter the discontinuity and
-        return oo:
+        return infinity:
 
         >>> e.as_sum(5, 'left')
-        oo
-        >>> e.as_sum(5, 'trapezoid')
-        oo
+        zoo
+
+        The number of intervals can be left unspecified.
+        >>> e = Integral(x**2, (x, 0, 2))
+        >>> e.as_sum(n, 'right').expand()
+        8/3 + 4/n + 4/(3*n**2)
+
+        This shows that the midpoint rule is more accurate, as its error
+        term decays as the square of n:
+        >>> e.as_sum(n, 'midpoint').expand()
+        8/3 - 2/(3*n**2)
+
+        A symbolic sum is returned with evaluate=False:
+        >>> e.as_sum(n, 'midpoint', evaluate=False)
+        2*Sum((2*_k/n - 1/n)**2, (_k, 1, n))/n
 
         See Also
         ========
@@ -1124,48 +1141,38 @@ def as_sum(self, n, method="midpoint"):
         Integral.doit : Perform the integration using any hints
         """
 
+        from sympy.concrete.summations import Sum
         limits = self.limits
         if len(limits) > 1:
             raise NotImplementedError(
                 "Multidimensional midpoint rule not implemented yet")
         else:
             limit = limits[0]
-            if len(limit) != 3:
-                raise ValueError("Expecting a definite integral.")
-        if n <= 0:
-            raise ValueError("n must be > 0")
-        if n == oo:
-            raise NotImplementedError("Infinite summation not yet implemented")
-        sym, lower_limit, upper_limit = limit
-        dx = (upper_limit - lower_limit)/n
-
-        if method == 'trapezoid':
-            l = self.function.limit(sym, lower_limit)
-            r = self.function.limit(sym, upper_limit, "-")
-            result = (l + r)/2
-            for i in range(1, n):
-                x = lower_limit + i*dx
-                result += self.function.subs(sym, x)
-            return result*dx
-        elif method not in ('left', 'right', 'midpoint'):
-            raise NotImplementedError("Unknown method %s" % method)
-
-        result = 0
-        for i in range(n):
-            if method == "midpoint":
-                xi = lower_limit + i*dx + dx/2
-            elif method == "left":
-                xi = lower_limit + i*dx
-                if i == 0:
-                    result = self.function.limit(sym, lower_limit)
-                    continue
-            elif method == "right":
-                xi = lower_limit + i*dx + dx
-                if i == n:
-                    result += self.function.limit(sym, upper_limit, "-")
-                    continue
-            result += self.function.subs(sym, xi)
-        return result*dx
+            if (len(limit) != 3 or limit[1].is_finite is False or
+                limit[2].is_finite is False):
+                raise ValueError("Expecting a definite integral over "
+                                  "a finite interval.")
+        n = sympify(n)
+        if (n.is_positive is False or n.is_integer is False or
+            n.is_finite is False):
+            raise ValueError("n must be a positive integer, got %s" % n)
+        x, a, b = limit
+        dx = (b - a)/n
+        k = Dummy('k', integer=True, positive=True)
+        f = self.function
+
+        if method == "left":
+            result = dx*Sum(f.subs(x, a + (k-1)*dx), (k, 1, n))
+        elif method == "right":
+            result = dx*Sum(f.subs(x, a + k*dx), (k, 1, n))
+        elif method == "midpoint":
+            result = dx*Sum(f.subs(x, a + k*dx - dx/2), (k, 1, n))
+        elif method == "trapezoid":
+            result = dx*((f.subs(x, a) + f.subs(x, b))/2 +
+                Sum(f.subs(x, a + k*dx), (k, 1, n - 1)))
+        else:
+            raise ValueError("Unknown method %s" % method)
+        return result.doit() if evaluate else result
 
     def _sage_(self):
         import sage.all as sage
diff --git a/sympy/integrals/tests/test_integrals.py b/sympy/integrals/tests/test_integrals.py
index bdab321a8caa..b13bea5ab5cd 100644
--- a/sympy/integrals/tests/test_integrals.py
+++ b/sympy/integrals/tests/test_integrals.py
@@ -683,10 +683,13 @@ def test_as_sum_midpoint1():
 
 def test_as_sum_midpoint2():
     e = Integral((x + y)**2, (x, 0, 1))
+    n = Symbol('n', positive=True, integer=True)
     assert e.as_sum(1, method="midpoint").expand() == S(1)/4 + y + y**2
     assert e.as_sum(2, method="midpoint").expand() == S(5)/16 + y + y**2
     assert e.as_sum(3, method="midpoint").expand() == S(35)/108 + y + y**2
     assert e.as_sum(4, method="midpoint").expand() == S(21)/64 + y + y**2
+    assert e.as_sum(n, method="midpoint").expand() == \
+        y**2 + y + 1/3 - 1/(12*n**2)
 
 
 def test_as_sum_left():
@@ -695,7 +698,9 @@ def test_as_sum_left():
     assert e.as_sum(2, method="left").expand() == S(1)/8 + y/2 + y**2
     assert e.as_sum(3, method="left").expand() == S(5)/27 + 2*y/3 + y**2
     assert e.as_sum(4, method="left").expand() == S(7)/32 + 3*y/4 + y**2
-
+    assert e.as_sum(n, method="left").expand() == \
+        y**2 + y + S(1)/3 - y/n - 1/(2*n) + 1/(6*n**2)
+    assert e.as_sum(10, method="left", evaluate=False).has(Sum)
 
 def test_as_sum_right():
     e = Integral((x + y)**2, (x, 0, 1))
@@ -703,15 +708,27 @@ def test_as_sum_right():
     assert e.as_sum(2, method="right").expand() == S(5)/8 + 3*y/2 + y**2
     assert e.as_sum(3, method="right").expand() == S(14)/27 + 4*y/3 + y**2
     assert e.as_sum(4, method="right").expand() == S(15)/32 + 5*y/4 + y**2
+    assert e.as_sum(n, method="right").expand() == \
+        y**2 + y + S(1)/3 + y/n + 1/(2*n) + 1/(6*n**2)
+
 
+def test_as_sum_trapezoid():
+    e = Integral((x + y)**2, (x, 0, 1))
+    assert e.as_sum(1, method="trapezoid").expand() == y**2 + y + S(1)/2
+    assert e.as_sum(2, method="trapezoid").expand() == y**2 + y + S(3)/8
+    assert e.as_sum(3, method="trapezoid").expand() == y**2 + y + S(19)/54
+    assert e.as_sum(4, method="trapezoid").expand() == y**2 + y + S(11)/32
+    assert e.as_sum(n, method="trapezoid").expand() == \
+        y**2 + y + S(1)/3 + 1/(6*n**2)
+    assert Integral(sign(x), (x, 0, 1)).as_sum(1, 'trapezoid') == S(1)/2
 
 def test_as_sum_raises():
     e = Integral((x + y)**2, (x, 0, 1))
     raises(ValueError, lambda: e.as_sum(-1))
     raises(ValueError, lambda: e.as_sum(0))
     raises(ValueError, lambda: Integral(x).as_sum(3))
-    raises(NotImplementedError, lambda: e.as_sum(oo))
-    raises(NotImplementedError, lambda: e.as_sum(3, method='xxxx2'))
+    raises(ValueError, lambda: e.as_sum(oo))
+    raises(ValueError, lambda: e.as_sum(3, method='xxxx2'))
 
 
 def test_nested_doit():


Code review comments:
diff_hunk: @@ -1047,23 +1047,29 @@ def _eval_as_leading_term(self, x):
                 break
         return integrate(leading_term, *self.args[1:])
 
-    def as_sum(self, n, method="midpoint"):
+    def as_sum(self, n, method="midpoint", evaluate=True):
path: sympy/integrals/integrals.py
text: @user1:
How about using `Symbol('n', integer=True, positive=True)` as default `n` if `n` is not given?

@author:
Yes, it makes sense. But I guess it should be a dummy if a symbol isn't provided, to avoid clashing with any `n` in user's context.