Repository: sympy/sympy
Pull request: Problem with solving simple separable ODE fixed
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->


#### Brief description of what is fixed or changed
Fixes #22155

#### Other comments


#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
NO ENTRY
<!-- END RELEASE NOTES -->


Problem statement:
Problem with solving simple separable ODE
I want to solve by SymPy simple separable equation y'=e^(y-x), but SymPy can't solve it!
```
from sympy import *
x = symbols("x")
y = Function("y")
eq = Eq(Derivative(y(x), x), exp(y(x)-x))
dsolve(eq, y(x))
```
But it solves successfuly similar equation y'=e^(y+x).

Follоwing code
```
eq = Eq(Derivative(y(x), x), exp(y(x)-x))
print(classify_ode(eq,y(x)))
```
returns

`('factorable', 'separable', '1st_power_series', 'lie_group', 'separable_Integral')
`

The method for 'factorable' seems broken. Explicitly using the second method,

`dsolve(eq, y(x), hint='separable')
`

quickly returns the correct solution.

For the second equation y'=e^(y+x), the list of methods is

`('separable', '1st_power_series', 'lie_group', 'separable_Integral')
`

so in this case the problematic method 'factorable' is not used. 

Experiments have shown that those equations cannot be solved if their exponents have a minus sign.

Hints from issue discussion:
The problem is that the factorable hint calls `dsolve` recursively but for this example ODE the recursive call will also end up being handled by factorable again leading to infinite recursion. It looks like the last line here is intended to prevent the infinite recursion:
https://github.com/sympy/sympy/blob/19b187982f90c868b987a8bbea085874647c98aa/sympy/solvers/ode/single.py#L897-L902
That check fails though because the expressions being compared are unequal:
```
903  ->	                return fraction(factor(self.eqs[0]))[0]-eq!=0
904  	            return True
905  	        for i in factors:
906  	            if i.has(f(x)):
907  	                self.eqs.append(i)
908  	        return len(self.eqs)>0 and len(factors)>1
(Pdb) p self.eqs[0]
-exp(-x + y(x)) + Derivative(y(x), x)
(Pdb) p fraction(factor(self.eqs[0]))[0]
Derivative(y(x), x) - exp(-x)*exp(y(x))
(Pdb) p eq
exp(x)*Derivative(y(x), x) - exp(y(x))
```
The ultimate fix is that factorable should not be a solver hint and should just be part of the preprocessing done by dsolve. Maybe for now that check could be improved. This diff fixes the example shown:
```diff
diff --git a/sympy/solvers/ode/single.py b/sympy/solvers/ode/single.py
index 7523b45..a3e5c42 100644
--- a/sympy/solvers/ode/single.py
+++ b/sympy/solvers/ode/single.py
@@ -892,14 +892,8 @@ def _matches(self):
             if len(self.eqs)>0:
                 return True
         roots = solve(eq, df)
-        if len(roots)>0:
+        if len(roots)>1:
             self.eqs = [(df - root) for root in roots]
-            if len(self.eqs)==1:
-                if order>1:
-                    return False
-                if self.eqs[0].has(Float):
-                    return False
-                return fraction(factor(self.eqs[0]))[0]-eq!=0
             return True
         for i in factors:
             if i.has(f(x)):

```
That might potentialy cause problems in other cases though. I'm not sure I understand the logic of this code.
By the way [online version](https://live.sympy.org/) of SymPy (version 1.5.1) solves this kind of equations correctly! 
I would like to help
> I would like to help

You could try applying the diff I showed as a PR to see if it breaks any of the tests.
> By the way [online version](https://live.sympy.org/) of SymPy (version 1.5.1) solves this kind of equations correctly!

Using `git bisect` shows that this was working until 1a26cdd3230b543775c77ba2493500166c0ff984 from  #18623
CC @Mohitbalwani26 
> CC @Mohitbalwani26

Looking into it. Thanks
@oscarbenjamin  Sure
@mohajain are you trying this?
Also i think this diff will work : 
```diff
--- a/sympy/solvers/ode/single.py
+++ b/sympy/solvers/ode/single.py
@@ -875,13 +875,13 @@ class Factorable(SingleODESolver):
     has_integral = False

     def _matches(self):
-        eq = self.ode_problem.eq
+        orig_eq = self.ode_problem.eq
         f = self.ode_problem.func.func
         x = self.ode_problem.sym
         order =self.ode_problem.order
         df = f(x).diff(x)
         self.eqs = []
-        eq = eq.collect(f(x), func = cancel)
+        eq = orig_eq.collect(f(x), func = cancel)
         eq = fraction(factor(eq))[0]
         factors = Mul.make_args(factor(eq))
         roots = [fac.as_base_exp() for fac in factors if len(fac.args)!=0]
@@ -899,7 +899,7 @@ def _matches(self):
                     return False
                 if self.eqs[0].has(Float):
                     return False
-                return fraction(factor(self.eqs[0]))[0]-eq!=0
+                return fraction(factor(self.eqs[0]))[0]-eq!=0 and simplify(orig_eq - fraction(factor(self.eqs[0]))[0])!=0
             return True
         for i in factors:
             if i.has(f(x)):
```
Using simplify like that seems flakey. There needs to be a clear robust way of preventing infinite recursion here that does not depend on simplifying anything.

Commit message: test case added

Patch to review:
diff --git a/sympy/solvers/ode/ode.py b/sympy/solvers/ode/ode.py
index 5bc8b7d65273..d9f55bf4acc0 100644
--- a/sympy/solvers/ode/ode.py
+++ b/sympy/solvers/ode/ode.py
@@ -931,7 +931,7 @@ class in it.  Note that a hint may do this anyway if
     '1st_homogeneous_coeff_subs_indep_div_dep_Integral',
     '1st_homogeneous_coeff_subs_dep_div_indep_Integral')
     >>> classify_ode(f(x).diff(x, 2) + 3*f(x).diff(x) + 2*f(x) - 4)
-    ('nth_linear_constant_coeff_undetermined_coefficients',
+    ('factorable', 'nth_linear_constant_coeff_undetermined_coefficients',
     'nth_linear_constant_coeff_variation_of_parameters',
     'nth_linear_constant_coeff_variation_of_parameters_Integral')
 
diff --git a/sympy/solvers/ode/single.py b/sympy/solvers/ode/single.py
index 7523b4536e86..1076c63ccc9b 100644
--- a/sympy/solvers/ode/single.py
+++ b/sympy/solvers/ode/single.py
@@ -12,7 +12,7 @@
 from sympy.core.exprtools import factor_terms
 from sympy.core.expr import Expr
 from sympy.core.function import AppliedUndef, Derivative, diff, Function, expand, Subs, _mexpand
-from sympy.core.numbers import Float, zoo
+from sympy.core.numbers import zoo
 from sympy.core.relational import Equality, Eq
 from sympy.core.symbol import Symbol, Dummy, Wild
 from sympy.core.mul import Mul
@@ -875,13 +875,12 @@ class Factorable(SingleODESolver):
     has_integral = False
 
     def _matches(self):
-        eq = self.ode_problem.eq
+        eq_orig = self.ode_problem.eq
         f = self.ode_problem.func.func
         x = self.ode_problem.sym
-        order =self.ode_problem.order
         df = f(x).diff(x)
         self.eqs = []
-        eq = eq.collect(f(x), func = cancel)
+        eq = eq_orig.collect(f(x), func = cancel)
         eq = fraction(factor(eq))[0]
         factors = Mul.make_args(factor(eq))
         roots = [fac.as_base_exp() for fac in factors if len(fac.args)!=0]
@@ -894,13 +893,9 @@ def _matches(self):
         roots = solve(eq, df)
         if len(roots)>0:
             self.eqs = [(df - root) for root in roots]
-            if len(self.eqs)==1:
-                if order>1:
-                    return False
-                if self.eqs[0].has(Float):
-                    return False
-                return fraction(factor(self.eqs[0]))[0]-eq!=0
-            return True
+            # Avoid infinite recursion
+            matches = self.eqs != [eq_orig]
+            return matches
         for i in factors:
             if i.has(f(x)):
                 self.eqs.append(i)
diff --git a/sympy/solvers/ode/tests/test_ode.py b/sympy/solvers/ode/tests/test_ode.py
index 1a92901bf8ca..a90dbf970e3c 100644
--- a/sympy/solvers/ode/tests/test_ode.py
+++ b/sympy/solvers/ode/tests/test_ode.py
@@ -96,21 +96,21 @@ def test_dsolve_options():
         '1st_homogeneous_coeff_subs_indep_div_dep_Integral', '1st_linear',
         '1st_linear_Integral', 'Bernoulli', 'Bernoulli_Integral',
         'almost_linear', 'almost_linear_Integral', 'best', 'best_hint',
-        'default', 'lie_group',
+        'default', 'factorable', 'lie_group',
         'nth_linear_euler_eq_homogeneous', 'order',
         'separable', 'separable_Integral']
     Integral_keys = ['1st_exact_Integral',
         '1st_homogeneous_coeff_subs_dep_div_indep_Integral',
         '1st_homogeneous_coeff_subs_indep_div_dep_Integral', '1st_linear_Integral',
         'Bernoulli_Integral', 'almost_linear_Integral', 'best', 'best_hint', 'default',
-        'nth_linear_euler_eq_homogeneous',
+        'factorable', 'nth_linear_euler_eq_homogeneous',
         'order', 'separable_Integral']
     assert sorted(a.keys()) == keys
     assert a['order'] == ode_order(eq, f(x))
     assert a['best'] == Eq(f(x), C1/x)
     assert dsolve(eq, hint='best') == Eq(f(x), C1/x)
-    assert a['default'] == 'separable'
-    assert a['best_hint'] == 'separable'
+    assert a['default'] == 'factorable'
+    assert a['best_hint'] == 'factorable'
     assert not a['1st_exact'].has(Integral)
     assert not a['separable'].has(Integral)
     assert not a['1st_homogeneous_coeff_best'].has(Integral)
@@ -126,8 +126,8 @@ def test_dsolve_options():
     assert b['order'] == ode_order(eq, f(x))
     assert b['best'] == Eq(f(x), C1/x)
     assert dsolve(eq, hint='best', simplify=False) == Eq(f(x), C1/x)
-    assert b['default'] == 'separable'
-    assert b['best_hint'] == '1st_linear'
+    assert b['default'] == 'factorable'
+    assert b['best_hint'] == 'factorable'
     assert a['separable'] != b['separable']
     assert a['1st_homogeneous_coeff_subs_dep_div_indep'] != \
         b['1st_homogeneous_coeff_subs_dep_div_indep']
@@ -231,7 +231,8 @@ def test_classify_ode():
          '1st_linear_Integral',
          'Bernoulli_Integral',
          'nth_linear_constant_coeff_variation_of_parameters_Integral')
-    assert c == ('1st_linear',
+    assert c == ('factorable',
+         '1st_linear',
          'Bernoulli',
          '1st_power_series',
          'lie_group',
@@ -243,7 +244,7 @@ def test_classify_ode():
 
     assert classify_ode(
         2*x*f(x)*f(x).diff(x) + (1 + x)*f(x)**2 - exp(x), f(x)
-    ) == ('1st_exact', 'Bernoulli', 'almost_linear', 'lie_group',
+    ) == ('factorable', '1st_exact', 'Bernoulli', 'almost_linear', 'lie_group',
         '1st_exact_Integral', 'Bernoulli_Integral', 'almost_linear_Integral')
     assert 'Riccati_special_minus2' in \
         classify_ode(2*f(x).diff(x) + f(x)**2 - f(x)/x + 3*x**(-2), f(x))
@@ -253,11 +254,11 @@ def test_classify_ode():
     k = Symbol('k')
     assert classify_ode(f(x).diff(x)/(k*f(x) + k*x*f(x)) + 2*f(x)/(k*f(x) +
         k*x*f(x)) + x*f(x).diff(x)/(k*f(x) + k*x*f(x)) + z, f(x)) == \
-        ('separable', '1st_exact', '1st_linear', 'Bernoulli',
+        ('factorable', 'separable', '1st_exact', '1st_linear', 'Bernoulli',
         '1st_power_series', 'lie_group', 'separable_Integral', '1st_exact_Integral',
         '1st_linear_Integral', 'Bernoulli_Integral')
     # preprocessing
-    ans = ('nth_algebraic', 'separable', '1st_exact', '1st_linear', 'Bernoulli',
+    ans = ('factorable', 'nth_algebraic', 'separable', '1st_exact', '1st_linear', 'Bernoulli',
         '1st_homogeneous_coeff_best',
         '1st_homogeneous_coeff_subs_indep_div_dep',
         '1st_homogeneous_coeff_subs_dep_div_indep',
@@ -724,7 +725,7 @@ def test_undetermined_coefficients_match():
 def test_issue_4785():
     from sympy.abc import A
     eq = x + A*(x + diff(f(x), x) + f(x)) + diff(f(x), x) + f(x) + 2
-    assert classify_ode(eq, f(x)) == ('1st_exact', '1st_linear',
+    assert classify_ode(eq, f(x)) == ('factorable', '1st_exact', '1st_linear',
         'almost_linear', '1st_power_series', 'lie_group',
         'nth_linear_constant_coeff_undetermined_coefficients',
         'nth_linear_constant_coeff_variation_of_parameters',
@@ -732,7 +733,7 @@ def test_issue_4785():
         'nth_linear_constant_coeff_variation_of_parameters_Integral')
     # issue 4864
     eq = (x**2 + f(x)**2)*f(x).diff(x) - 2*x*f(x)
-    assert classify_ode(eq, f(x)) == ('1st_exact',
+    assert classify_ode(eq, f(x)) == ('factorable', '1st_exact',
         '1st_homogeneous_coeff_best',
         '1st_homogeneous_coeff_subs_indep_div_dep',
         '1st_homogeneous_coeff_subs_dep_div_indep',
@@ -895,7 +896,7 @@ def test_2nd_power_series_ordinary():
     assert checkodesol(eq, sol) == (True, 0)
 
     eq = (1 + x**2)*(f(x).diff(x, 2)) + 2*x*(f(x).diff(x)) -2*f(x)
-    assert classify_ode(eq) == ('2nd_hypergeometric', '2nd_hypergeometric_Integral',
+    assert classify_ode(eq) == ('factorable', '2nd_hypergeometric', '2nd_hypergeometric_Integral',
     '2nd_power_series_ordinary')
 
     sol = Eq(f(x), C2*(-x**4/3 + x**2 + 1) + C1*x + O(x**6))
@@ -903,7 +904,7 @@ def test_2nd_power_series_ordinary():
     assert checkodesol(eq, sol) == (True, 0)
 
     eq = f(x).diff(x, 2) + x*(f(x).diff(x)) + f(x)
-    assert classify_ode(eq) == ('2nd_power_series_ordinary',)
+    assert classify_ode(eq) == ('factorable', '2nd_power_series_ordinary',)
     sol = Eq(f(x), C2*(x**4/8 - x**2/2 + 1) + C1*x*(-x**2/3 + 1) + O(x**6))
     assert dsolve(eq) == sol
     # FIXME: checkodesol fails for this solution...
diff --git a/sympy/solvers/ode/tests/test_single.py b/sympy/solvers/ode/tests/test_single.py
index 9b029d23d48c..ce422ba68d80 100644
--- a/sympy/solvers/ode/tests/test_single.py
+++ b/sympy/solvers/ode/tests/test_single.py
@@ -389,7 +389,7 @@ def test_lie_group():
 def test_separable_reduced():
     df = f(x).diff(x)
     eq = (x / f(x))*df  + tan(x**2*f(x) / (x**2*f(x) - 1))
-    assert classify_ode(eq) == ('separable_reduced', 'lie_group',
+    assert classify_ode(eq) == ('factorable', 'separable_reduced', 'lie_group',
         'separable_reduced_Integral')
     _ode_solver_test(_get_examples_ode_sol_separable_reduced)
 
@@ -1812,6 +1812,12 @@ def _get_examples_ode_sol_separable():
         'sol': [Eq(v(t), -68.585712797928991/tanh(C1 - 0.14288690166235204*t))],
         'func': v(t),
         'checkodesol_XFAIL': True,
+    },
+
+    #https://github.com/sympy/sympy/issues/22155
+    'separable_27': {
+        'eq': f(x).diff(x) - exp(f(x) - x),
+        'sol': [Eq(f(x), log(-exp(x)/(C1*exp(x) - 1)))],
     }
     }
     }
diff --git a/sympy/solvers/tests/test_pde.py b/sympy/solvers/tests/test_pde.py
index 1b43eb0b0886..da9e09ba18f6 100644
--- a/sympy/solvers/tests/test_pde.py
+++ b/sympy/solvers/tests/test_pde.py
@@ -231,4 +231,4 @@ def test_pdsolve_variable_coeff():
 
     eq = exp(2*x)*(u.diff(y)) + y*u - u
     sol = pdsolve(eq, hint='1st_linear_variable_coeff')
-    assert sol == Eq(u, exp((-y**2 + 2*y + 2*F(x))*exp(-2*x)/2))
+    assert sol == Eq(f(x, y), F(x)*exp(-y*(y - 2)*exp(-2*x)/2))


Code review comments:
diff_hunk: @@ -231,4 +231,4 @@ def test_pdsolve_variable_coeff():
 
     eq = exp(2*x)*(u.diff(y)) + y*u - u
     sol = pdsolve(eq, hint='1st_linear_variable_coeff')
-    assert sol == Eq(u, exp((-y**2 + 2*y + 2*F(x))*exp(-2*x)/2))
+    assert sol == Eq(f(x, y), F(x)*exp(-y*(y - 2)*exp(-2*x)/2))
path: sympy/solvers/tests/test_pde.py
text: @user1:
This should still be `u` rather than `f(x, y)` to keep it consistent with the examples above