Repository: sympy/sympy
Pull request: Rescale ODEs to match Euler solver
<!-- 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://github.com/blog/1506-closing-issues-via-pull-requests .-->

Fixes #15217 

Needs a fix for #15433 for one of the tests

#### Brief description of what is fixed or changed

When matching an ODE against the Euler form this patch multiplies the coefficients by an appropriate power of `x`. This means that `dsolve` can recognise more ODEs as being Euler ODEs and can solve some ODEs that it couldn't before.

This example gives `NotImplementedError` on current master:
```julia
In [7]: eqn = -exp(x) + (x*Derivative(f(x), (x, 2)) + Derivative(f(x), x))/x

In [8]: eqn
Out[8]: 
           2                 
          d          d       
       x⋅───(f(x)) + ──(f(x))
           2         dx      
   x     dx                  
- ℯ  + ──────────────────────
                 x           

In [9]: dsolve(eqn, f(x))
Out[9]: 
                                            ⌠               
                                 x          ⎮    x          
f(x) = C₁ + C₂⋅log(x) + (x - 1)⋅ℯ ⋅log(x) - ⎮ x⋅ℯ ⋅log(x) dx
                                            ⌡    
```

#### Other comments

I've had to adjust the tests in `classify_ode` because the Euler solver now matches many trivial ODEs.

I found #15433 while working on this. One of the tests uses an awkward workaround for that. Probably #15433 should be fixed first.

#### Release Notes

<!-- Write the release notes for this release below. 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 -->
* solvers
    * dsolve now rescales ODEs so that the Euler solver can be used in more cases
<!-- END RELEASE NOTES -->


Problem statement:
Euler ODE solver should recognize when there is a missing factor of x
```
>>> classify_ode(eqn.doit().expand())
()
>>> classify_ode((x**2*eqn.doit()).expand())
('nth_linear_euler_eq_nonhomogeneous_variation_of_parameters', 'nth_linear_euler_eq_nonhomogeneous_variation_of_parameters_Integral')
```

Probably the simplest way to fix this would be to make the ODE monic in the highest derivative and change the matching to check against that.

Hints from issue discussion:
What is `eqn` defined as?
Oh derp. It's from the [mailing list](https://groups.google.com/forum/#!topic/sympy/VBRmeCCAN0U). 

```
>>> eqn.doit()
-exp(x) + (x*Derivative(f(x), (x, 2)) + Derivative(f(x), x))/x
```

The original `eqn` is in factored form, but that's a separate issue from this (although I do think classify_ode should try expanding factored expressions and seeing if they can be classified, see also https://github.com/sympy/sympy/issues/6259).

```
>>> eqn
-exp(x) + Derivative(x*Derivative(f(x), x), x)/x
```

Commit message: Update classify_ode doctest after nth_euler fixes

Patch to review:
diff --git a/sympy/solvers/ode.py b/sympy/solvers/ode.py
index d95efe24498a..6f0a5e06bfdc 100644
--- a/sympy/solvers/ode.py
+++ b/sympy/solvers/ode.py
@@ -934,7 +934,8 @@ class in it.  Note that a hint may do this anyway if
     '1st_homogeneous_coeff_subs_indep_div_dep',
     '1st_homogeneous_coeff_subs_dep_div_indep',
     '1st_power_series', 'lie_group',
-    'nth_linear_constant_coeff_homogeneous', 'nth_algebraic_Integral',
+    'nth_linear_constant_coeff_homogeneous',
+    'nth_linear_euler_eq_homogeneous', 'nth_algebraic_Integral',
     'separable_Integral', '1st_linear_Integral',
     '1st_homogeneous_coeff_subs_indep_div_dep_Integral',
     '1st_homogeneous_coeff_subs_dep_div_indep_Integral')
@@ -1409,17 +1410,28 @@ def _test_term(coeff, order):
             elif order == 1:
                 return x == coeff
             return False
-        if r and not any(not _test_term(r[i], i) for i in r if i >= 0):
-            if not r[-1]:
-                matching_hints["nth_linear_euler_eq_homogeneous"] = r
+
+        # Find coefficient for higest derivative, multiply coefficients to
+        # bring the equation into Euler form if possible
+        r_rescaled = None
+        if r is not None:
+            coeff = r[order]
+            power = coeff.match(k*x**n)[n]
+            factor = x**(order - power)
+            r_rescaled = {i: factor*r[i] for i in r}
+
+        if r_rescaled and not any(not _test_term(r_rescaled[i], i) for i in
+                r_rescaled if i != 'trialset' and i >= 0):
+            if not r_rescaled[-1]:
+                matching_hints["nth_linear_euler_eq_homogeneous"] = r_rescaled
             else:
-                matching_hints["nth_linear_euler_eq_nonhomogeneous_variation_of_parameters"] = r
-                matching_hints["nth_linear_euler_eq_nonhomogeneous_variation_of_parameters_Integral"] = r
-                e, re = posify(r[-1].subs(x, exp(x)))
+                matching_hints["nth_linear_euler_eq_nonhomogeneous_variation_of_parameters"] = r_rescaled
+                matching_hints["nth_linear_euler_eq_nonhomogeneous_variation_of_parameters_Integral"] = r_rescaled
+                e, re = posify(r_rescaled[-1].subs(x, exp(x)))
                 undetcoeff = _undetermined_coefficients_match(e.subs(re), x)
                 if undetcoeff['test']:
-                    r['trialset'] = undetcoeff['trialset']
-                    matching_hints["nth_linear_euler_eq_nonhomogeneous_undetermined_coefficients"] = r
+                    r_rescaled['trialset'] = undetcoeff['trialset']
+                    matching_hints["nth_linear_euler_eq_nonhomogeneous_undetermined_coefficients"] = r_rescaled
 
 
     # Order keys based on allhints.
diff --git a/sympy/solvers/tests/test_ode.py b/sympy/solvers/tests/test_ode.py
index 76ca1fcfae7b..31ce2dff5ec9 100644
--- a/sympy/solvers/tests/test_ode.py
+++ b/sympy/solvers/tests/test_ode.py
@@ -717,6 +717,7 @@ def test_classify_ode():
     assert classify_ode(f(x).diff(x, 2), f(x)) == \
         ('nth_algebraic',
         'nth_linear_constant_coeff_homogeneous',
+        'nth_linear_euler_eq_homogeneous',
         'Liouville',
         '2nd_power_series_ordinary',
         'nth_algebraic_Integral',
@@ -731,6 +732,7 @@ def test_classify_ode():
         '1st_homogeneous_coeff_subs_dep_div_indep',
         '1st_power_series', 'lie_group',
         'nth_linear_constant_coeff_homogeneous',
+        'nth_linear_euler_eq_homogeneous',
         'nth_algebraic_Integral',
         'separable_Integral',
         '1st_linear_Integral',
@@ -776,14 +778,17 @@ def test_classify_ode():
         '1st_homogeneous_coeff_subs_dep_div_indep',
         '1st_power_series', 'lie_group',
         'nth_linear_constant_coeff_undetermined_coefficients',
+        'nth_linear_euler_eq_nonhomogeneous_undetermined_coefficients',
         'nth_linear_constant_coeff_variation_of_parameters',
+        'nth_linear_euler_eq_nonhomogeneous_variation_of_parameters',
         'nth_algebraic_Integral',
         'separable_Integral', '1st_exact_Integral',
         '1st_linear_Integral',
         'Bernoulli_Integral',
         '1st_homogeneous_coeff_subs_indep_div_dep_Integral',
         '1st_homogeneous_coeff_subs_dep_div_indep_Integral',
-        'nth_linear_constant_coeff_variation_of_parameters_Integral')
+        'nth_linear_constant_coeff_variation_of_parameters_Integral',
+        'nth_linear_euler_eq_nonhomogeneous_variation_of_parameters_Integral')
     #     w/o f(x) given
     assert classify_ode(diff(f(x) + x, x) + diff(f(x), x)) == ans
     #     w/ f(x) and prep=True
@@ -791,8 +796,11 @@ def test_classify_ode():
                         prep=True) == ans
 
     assert classify_ode(Eq(2*x**3*f(x).diff(x), 0), f(x)) == \
-        ('nth_algebraic', 'separable', '1st_linear', '1st_power_series', 'lie_group',
-         'nth_algebraic_Integral', 'separable_Integral', '1st_linear_Integral')
+        ('nth_algebraic', 'separable', '1st_linear', '1st_power_series',
+         'lie_group', 'nth_linear_euler_eq_homogeneous',
+         'nth_algebraic_Integral', 'separable_Integral',
+         '1st_linear_Integral')
+
     assert classify_ode(Eq(2*f(x)**3*f(x).diff(x), 0), f(x)) == \
         ('nth_algebraic', 'separable', '1st_power_series', 'lie_group',
                 'nth_algebraic_Integral', 'separable_Integral')
@@ -2420,6 +2428,16 @@ def test_nth_order_linear_euler_eq_nonhomogeneous_variation_of_parameters():
     assert dsolve(eq, f(x), hint=our_hint).rhs in (sol, sols)
     assert checkodesol(eq, sol, order=2, solve_for_func=False)[0]
 
+    eq = -exp(x) + (x*Derivative(f(x), (x, 2)) + Derivative(f(x), x))/x
+    sol = Eq(f(x), C1 + C2*log(x) + (x - 1)*exp(x)*log(x) - Integral(x*exp(x)*log(x), x))
+    assert our_hint in classify_ode(eq)
+    dsolve_sol = dsolve(eq, f(x), hint=our_hint)
+    from sympy.integrals.risch import NonElementaryIntegral
+    from sympy.integrals.integrals import Integral as _Integral
+    dsolve_sol = dsolve_sol.replace(NonElementaryIntegral, lambda *args: _Integral(*args))
+    assert dsolve_sol == sol
+    assert checkodesol(eq, sol, order=2, solve_for_func=False)[0]
+
 
 def test_issue_5095():
     f = Function('f')


Code review comments:
diff_hunk: @@ -2420,6 +2428,16 @@ def test_nth_order_linear_euler_eq_nonhomogeneous_variation_of_parameters():
     assert dsolve(eq, f(x), hint=our_hint).rhs in (sol, sols)
     assert checkodesol(eq, sol, order=2, solve_for_func=False)[0]
 
+    eq = -exp(x) + (x*Derivative(f(x), (x, 2)) + Derivative(f(x), x))/x
+    sol = Eq(f(x), C1 + C2*log(x) + (x - 1)*exp(x)*log(x) - Integral(x*exp(x)*log(x), x))
path: sympy/solvers/tests/test_ode.py
text: @user1:
The following will work without calling `.replace()`:
```
    eq = -exp(x) + (x*Derivative(f(x), (x, 2)) + Derivative(f(x), x))/x
    from sympy.integrals.risch import NonElementaryIntegral
    sol = Eq(f(x), C1 + C2*log(x) + (x - 1)*exp(x)*log(x) - NonElementaryIntegral(x*exp(x)*log(x), x))
    assert our_hint in classify_ode(eq)
    assert dsolve(eq, f(x), hint=our_hint) == sol
    assert checkodesol(eq, sol, order=2, solve_for_func=False)[0]
```