Repository: sympy/sympy

Issue:
differentiate_finite gives wrong results for expressions with embedded derivatives
The `differentiate_finite` function works great for expressions that do not include any derivatives (like expressions in the [Sympy Tutorial](https://docs.sympy.org/latest/tutorial/calculus.html)), however when more complex expressions are used the `differentiate_finite` discretizes expressions in a wrong way.

If we want to use `differentiate_finite` to make finite difference of expressions that including functions and derivatives, e.g., `expr = nu(y)*diff(U(y), y)`, one would write

```
from sympy import *
y = Symbol('y') 
dy, U, nu = symbols('dy U nu', cls=Function)
df_expr = differentiate_finite(nu(y)*diff(U(y), y), y, points=dy(y))
```

Note that I am also trying to make dy = dy(y), which is common when making finite differences on stretched grids. Executing the code above does not lead to proper discretization. The proper discratization for the above expression should be

```assert (df_true - df_expr).simplify() 
df_true = nu(y+dy(y)/2)*(U(y+dy(y))/dy(y+dy(y)/2)-U(y)/dy(y+dy(y)/2))/dy(y) - nu(y-dy(y)/2)*(U(y)/dy(y-dy(y)/2)-U(y-dy(y))/dy(y-dy(y)/2))/dy(y)
```

I am extremely sorry for the poor readability of these expressions. Nevertheless, if this was made a test `assert (df_expr - df_true).simplify() == 0` would fail.

The problem as I see, is that `differentiate_finite` replaces the embedded derivatives with finite derivatives at the wrong `x0` locations. In order to have correct discretization of expressions with embedded derivatives the `x0` for embedded derivatives should depend on the `x0` and `points` arguments of the outer derivative.

What do you think about it?

Discussion:

Patch:
diff --git a/doc/src/tutorial/calculus.rst b/doc/src/tutorial/calculus.rst
index f4163b04ed..5b2da66216 100644
--- a/doc/src/tutorial/calculus.rst
+++ b/doc/src/tutorial/calculus.rst
@@ -341,14 +341,6 @@ the ``differentiate_finite`` function:
     >>> differentiate_finite(f(x)*g(x))
     -f(x - 1/2)⋅g(x - 1/2) + f(x + 1/2)⋅g(x + 1/2)
 
-If we want to expand the intermediate derivative we may pass the
-flag ``evaluate=True``:
-
-    >>> differentiate_finite(f(x)*g(x), evaluate=True)
-    (-f(x - 1/2) + f(x + 1/2))⋅g(x) + (-g(x - 1/2) + g(x + 1/2))⋅f(x)
-
-This form however does not respect the product rule.
-
 If you already have a ``Derivative`` instance, you can use the
 ``as_finite_difference`` method to generate approximations of the
 derivative to arbitrary order:
diff --git a/sympy/calculus/finite_diff.py b/sympy/calculus/finite_diff.py
index ab15e039c6..f09fca7060 100644
--- a/sympy/calculus/finite_diff.py
+++ b/sympy/calculus/finite_diff.py
@@ -18,8 +18,11 @@
 """
 
 from sympy import Derivative, S
+from sympy.core.basic import preorder_traversal
 from sympy.core.compatibility import iterable
 from sympy.core.decorators import deprecated
+from sympy.core.function import Subs
+from sympy.utilities.exceptions import SymPyDeprecationWarning
 
 
 def finite_diff_weights(order, x_list, x0=S.One):
@@ -29,7 +32,6 @@ def finite_diff_weights(order, x_list, x0=S.One):
     0, 1, ..., up to ``order`` using a recursive formula. Order of accuracy
     is at least ``len(x_list) - order``, if ``x_list`` is defined correctly.
 
-
     Parameters
     ==========
 
@@ -123,7 +125,6 @@ def finite_diff_weights(order, x_list, x0=S.One):
     Note that, unless you plan on using approximations based on subsets of
     ``x_list``, the order of gridpoints does not matter.
 
-
     The capability to generate weights at arbitrary points can be
     used e.g. to minimize Runge's phenomenon by using Chebyshev nodes:
 
@@ -157,7 +158,6 @@ def finite_diff_weights(order, x_list, x0=S.One):
 
     sympy.calculus.finite_diff.apply_finite_diff
 
-
     References
     ==========
 
@@ -243,9 +243,9 @@ def apply_finite_diff(order, x_list, y_list, x0=S.Zero):
     >>> i = Idx('i')
     >>> x_list, y_list = zip(*[(x[i+j], y[i+j]) for j in range(-1,2)])
     >>> apply_finite_diff(1, x_list, y_list, x[i])
-    ((x[i + 1] - x[i])/(-x[i - 1] + x[i]) - 1)*y[i]/(x[i + 1] - x[i]) - \
-(x[i + 1] - x[i])*y[i - 1]/((x[i + 1] - x[i - 1])*(-x[i - 1] + x[i])) + \
-(-x[i - 1] + x[i])*y[i + 1]/((x[i + 1] - x[i - 1])*(x[i + 1] - x[i]))
+    ((x[i + 1] - x[i])/(-x[i - 1] + x[i]) - 1)*y[i]/(x[i + 1] - x[i]) -
+    (x[i + 1] - x[i])*y[i - 1]/((x[i + 1] - x[i - 1])*(-x[i - 1] + x[i])) +
+    (-x[i - 1] + x[i])*y[i + 1]/((x[i + 1] - x[i - 1])*(x[i + 1] - x[i]))
 
     Notes
     =====
@@ -315,7 +315,6 @@ def _as_finite_diff(derivative, points=1, x0=None, wrt=None):
         derivative is to be approximated for. If not provided it
         is required that the Derivative is ordinary. Default: ``None``.
 
-
     Examples
     ========
 
@@ -347,10 +346,9 @@ def _as_finite_diff(derivative, points=1, x0=None, wrt=None):
     >>> e, sq2 = exp(1), sqrt(2)
     >>> xl = [x-h, x+h, x+e*h]
     >>> as_finite_diff(f(x).diff(x, 1), xl, x+h*sq2)
-    2*h*((h + sqrt(2)*h)/(2*h) - (-sqrt(2)*h + h)/(2*h))*f(E*h + x)/\
-((-h + E*h)*(h + E*h)) + (-(-sqrt(2)*h + h)/(2*h) - \
-(-sqrt(2)*h + E*h)/(2*h))*f(-h + x)/(h + E*h) + \
-(-(h + sqrt(2)*h)/(2*h) + (-sqrt(2)*h + E*h)/(2*h))*f(h + x)/(-h + E*h)
+    2*h*((h + sqrt(2)*h)/(2*h) - (-sqrt(2)*h + h)/(2*h))*f(E*h + x)/((-h + E*h)*(h + E*h)) +
+    (-(-sqrt(2)*h + h)/(2*h) - (-sqrt(2)*h + E*h)/(2*h))*f(-h + x)/(h + E*h) +
+    (-(h + sqrt(2)*h)/(2*h) + (-sqrt(2)*h + E*h)/(2*h))*f(h + x)/(-h + E*h)
 
     Partial derivatives are also supported:
 
@@ -430,18 +428,15 @@ def differentiate_finite(expr, *symbols,
 
     Parameters
     ==========
+
     expr : expression
     \*symbols : differentiate with respect to symbols
-    points: sequence or coefficient, optional
+    points: sequence, coefficient or undefined function, optional
         see ``Derivative.as_finite_difference``
     x0: number or Symbol, optional
         see ``Derivative.as_finite_difference``
     wrt: Symbol, optional
         see ``Derivative.as_finite_difference``
-    evaluate : bool
-        kwarg passed on to ``diff``, whether or not to
-        evaluate the Derivative intermediately (default: ``False``).
-
 
     Examples
     ========
@@ -452,21 +447,25 @@ def differentiate_finite(expr, *symbols,
     >>> differentiate_finite(f(x)*g(x), x, points=[x-h, x+h])
     -f(-h + x)*g(-h + x)/(2*h) + f(h + x)*g(h + x)/(2*h)
 
-    Note that the above form preserves the product rule in discrete form.
-    If we want we can pass ``evaluate=True`` to get another form (which is
-    usually not what we want):
-
-    >>> differentiate_finite(f(x)*g(x), x, points=[x-h, x+h], evaluate=True).simplify()
-    -((f(-h + x) - f(h + x))*g(x) + (g(-h + x) - g(h + x))*f(x))/(2*h)
-
-    ``differentiate_finite`` works on any expression:
+    ``differentiate_finite`` works on any expression, including the expressions
+    with embedded derivatives:
 
     >>> differentiate_finite(f(x) + sin(x), x, 2)
     -2*f(x) + f(x - 1) + f(x + 1) - 2*sin(x) + sin(x - 1) + sin(x + 1)
-    >>> differentiate_finite(f(x) + sin(x), x, 2, evaluate=True)
-    -2*f(x) + f(x - 1) + f(x + 1) - sin(x)
     >>> differentiate_finite(f(x, y), x, y)
     f(x - 1/2, y - 1/2) - f(x - 1/2, y + 1/2) - f(x + 1/2, y - 1/2) + f(x + 1/2, y + 1/2)
+    >>> differentiate_finite(f(x)*g(x).diff(x), x)
+    (-g(x) + g(x + 1))*f(x + 1/2) - (g(x) - g(x - 1))*f(x - 1/2)
+
+    To make finite difference with non-constant discretization step use
+    undefined functions:
+
+    >>> dx = Function('dx')
+    >>> differentiate_finite(f(x)*g(x).diff(x), points=dx(x))
+    -(-g(x - dx(x)/2 - dx(x - dx(x)/2)/2)/dx(x - dx(x)/2) +
+    g(x - dx(x)/2 + dx(x - dx(x)/2)/2)/dx(x - dx(x)/2))*f(x - dx(x)/2)/dx(x) +
+    (-g(x + dx(x)/2 - dx(x + dx(x)/2)/2)/dx(x + dx(x)/2) +
+    g(x + dx(x)/2 + dx(x + dx(x)/2)/2)/dx(x + dx(x)/2))*f(x + dx(x)/2)/dx(x)
 
     """
     # Key-word only arguments only available in Python 3
@@ -474,10 +473,22 @@ def differentiate_finite(expr, *symbols,
     x0 = kwargs.pop('x0', None)
     wrt = kwargs.pop('wrt', None)
     evaluate = kwargs.pop('evaluate', False)
+    if any(term.is_Derivative for term in list(preorder_traversal(expr))):
+        evaluate = False
     if kwargs:
         raise ValueError("Unknown kwargs: %s" % kwargs)
 
     Dexpr = expr.diff(*symbols, evaluate=evaluate)
-    return Dexpr.replace(
-        lambda arg: arg.is_Derivative,
-        lambda arg: arg.as_finite_difference(points=points, x0=x0, wrt=wrt))
+    if evaluate:
+        SymPyDeprecationWarning(feature="``evaluate`` flag",
+                                issue=17881,
+                                deprecated_since_version="1.5").warn()
+        return Dexpr.replace(
+            lambda arg: arg.is_Derivative,
+            lambda arg: arg.as_finite_difference(points=points, x0=x0, wrt=wrt))
+    else:
+        DFexpr = Dexpr.as_finite_difference(points=points, x0=x0, wrt=wrt)
+        return DFexpr.replace(
+            lambda arg: isinstance(arg, Subs),
+            lambda arg: arg.expr.as_finite_difference(
+                    points=points, x0=arg.point[0], wrt=arg.variables[0]))


Test patch:
diff --git a/sympy/calculus/tests/test_finite_diff.py b/sympy/calculus/tests/test_finite_diff.py
index 868fac0587..4a53a3a33d 100644
--- a/sympy/calculus/tests/test_finite_diff.py
+++ b/sympy/calculus/tests/test_finite_diff.py
@@ -5,7 +5,8 @@
     apply_finite_diff, differentiate_finite, finite_diff_weights,
     as_finite_diff
 )
-from sympy.testing.pytest import raises, warns_deprecated_sympy
+from sympy.testing.pytest import raises, warns_deprecated_sympy, ignore_warnings
+from sympy.utilities.exceptions import SymPyDeprecationWarning
 
 
 def test_apply_finite_diff():
@@ -123,15 +124,17 @@ def test_as_finite_diff():
 
 
 def test_differentiate_finite():
-    x, y = symbols('x y')
+    x, y, h = symbols('x y h')
     f = Function('f')
-    res0 = differentiate_finite(f(x, y) + exp(42), x, y, evaluate=True)
+    with ignore_warnings(SymPyDeprecationWarning):
+        res0 = differentiate_finite(f(x, y) + exp(42), x, y, evaluate=True)
     xm, xp, ym, yp = [v + sign*S.Half for v, sign in product([x, y], [-1, 1])]
     ref0 = f(xm, ym) + f(xp, yp) - f(xm, yp) - f(xp, ym)
     assert (res0 - ref0).simplify() == 0
 
     g = Function('g')
-    res1 = differentiate_finite(f(x)*g(x) + 42, x, evaluate=True)
+    with ignore_warnings(SymPyDeprecationWarning):
+        res1 = differentiate_finite(f(x)*g(x) + 42, x, evaluate=True)
     ref1 = (-f(x - S.Half) + f(x + S.Half))*g(x) + \
            (-g(x - S.Half) + g(x + S.Half))*f(x)
     assert (res1 - ref1).simplify() == 0
@@ -141,3 +144,23 @@ def test_differentiate_finite():
     assert (res2 - ref2).simplify() == 0
     raises(ValueError, lambda: differentiate_finite(f(x)*g(x), x,
                                                     pints=[x-1, x+1]))
+
+    res3 = differentiate_finite(f(x)*g(x).diff(x), x)
+    ref3 = (-g(x) + g(x + 1))*f(x + S.Half) - (g(x) - g(x - 1))*f(x - S.Half)
+    assert res3 == ref3
+
+    res4 = differentiate_finite(f(x)*g(x).diff(x).diff(x), x)
+    ref4 = -((g(x - Rational(3, 2)) - 2*g(x - S.Half) + g(x + S.Half))*f(x - S.Half)) \
+           + (g(x - S.Half) - 2*g(x + S.Half) + g(x + Rational(3, 2)))*f(x + S.Half)
+    assert res4 == ref4
+
+    res5_expr = f(x).diff(x)*g(x).diff(x)
+    res5 = differentiate_finite(res5_expr, points=[x-h, x, x+h])
+    ref5 = (-2*f(x)/h + f(-h + x)/(2*h) + 3*f(h + x)/(2*h))*(-2*g(x)/h + g(-h + x)/(2*h) \
+           + 3*g(h + x)/(2*h))/(2*h) - (2*f(x)/h - 3*f(-h + x)/(2*h) - \
+           f(h + x)/(2*h))*(2*g(x)/h - 3*g(-h + x)/(2*h) - g(h + x)/(2*h))/(2*h)
+    assert res5 == ref5
+
+    res6 = res5.limit(h, 0).doit()
+    ref6 = diff(res5_expr, x)
+    assert res6 == ref6