Repository: sympy/sympy
Pull request: dsolve gave NotImplemented error while solving with ics
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->
Fixes #15724
#### 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 . Please also
write a comment on that issue linking back to this pull request once it is
open. --> 
The main PR for the issue was #16098, but due to some issues It was closed and so now this PR contains all the changes that were meant to be done.

Earlier the output for dsolve (given ics) was

```
>>> dsolve(y(x).diff(x) - 1/y(x), y(x), ics={y(1): 2})
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/k3/sympy/sympy/solvers/ode.py", line 680, in dsolve
    return _helper_simplify(eq, hint, hints, simplify, ics=ics)
  File "/home/k3/sympy/sympy/solvers/ode.py", line 723, in _helper_simplify
    solved_constants = solve_ics([s], [r['func']], cons(s), ics)
  File "/home/k3/sympy/sympy/solvers/ode.py", line 818, in solve_ics
    raise NotImplementedError("Couldn't solve for initial conditions")
NotImplementedError: Couldn't solve for initial conditions
```

Now the output is,

```
>>> dsolve(y(x).diff(x) - 1/y(x), y(x), ics={y(1): 2})
[Eq(y(x), sqrt(2*x + 2))]
```

#### Brief description of what is fixed or changed


#### Other comments


#### 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 can now choose from multiple general solutions based in initial conditions.
<!-- END RELEASE NOTES -->


Problem statement:
dsolve fails to solve for initial condition which does not seem complicated
Trying to solve `y' = 1/y` with initial condition `y(1) = 2`. Expected solution: `y = sqrt(2*x + 2)`. Works fine without ICS, fails with it. 
```
>>> from sympy import *
>>> x, y = Symbol('x'), Function('y')
>>> dsolve(y(x).diff(x) - 1/y(x), y(x))
[Eq(y(x), -sqrt(C1 + 2*x)), Eq(y(x), sqrt(C1 + 2*x))]
>>> dsolve(y(x).diff(x) - 1/y(x), y(x), ics={y(1): 2})
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/k3/sympy/sympy/solvers/ode.py", line 680, in dsolve
    return _helper_simplify(eq, hint, hints, simplify, ics=ics)
  File "/home/k3/sympy/sympy/solvers/ode.py", line 723, in _helper_simplify
    solved_constants = solve_ics([s], [r['func']], cons(s), ics)
  File "/home/k3/sympy/sympy/solvers/ode.py", line 818, in solve_ics
    raise NotImplementedError("Couldn't solve for initial conditions")
NotImplementedError: Couldn't solve for initial conditions
```

Hints from issue discussion:
There could be a try / catch block on line 723 in `ode.py`. `solve_ics` can raise 4 different errors (either ValueError or NotImplementedError) that one might like to do different things with- or just ignore all of them and `pass`. Right now if any general solution cannot be solved (for initial conditions) it won't return any of them.
The IC code isn't properly setup to handle the case of multiple general solutions. In this case one form works and the other does not:
```julia
In [12]: solve(sqrt(x+2)-2)                                                                                                                                   
Out[12]: [2]

In [13]: solve(-sqrt(x+2)-2)                                                                                                                                  
Out[13]: []
```
In the second case a simpler example is
```julia
In [18]: solve(Eq(sqrt(x), -2))                                                                                                                               
Out[18]: []
```
The code for handling ICs should consider that when there are multiple solution forms the ICs can be used to select the valid form. In this case one form is valid for positive ICs and one for negative ICs.
@oscarbenjamin @ethankward Can you please review the PR ?
I think this might be the place that normally handles that but it should really be the very last thing that happens in `dsolve`:
[https://github.com/sympy/sympy/blob/fe3dea79551bbeab86011361af2364abc53105f6/sympy/solvers/ode.py#L2330](https://github.com/sympy/sympy/blob/fe3dea79551bbeab86011361af2364abc53105f6/sympy/solvers/ode.py#L2330)
Sorry that last message was intended for the PR. I've recommented there.

Commit message: Corrected assertion error

Patch to review:
diff --git a/sympy/solvers/ode.py b/sympy/solvers/ode.py
index 2b66aa1baee5..5998ca5bfc08 100644
--- a/sympy/solvers/ode.py
+++ b/sympy/solvers/ode.py
@@ -714,8 +714,13 @@ def _helper_simplify(eq, hint, match, simplify=True, ics=None, **kwargs):
         else:
             rv1 = []
             for s in rv:
-                solved_constants = solve_ics([s], [r['func']], cons(s), ics)
+                try:
+                    solved_constants = solve_ics([s], [r['func']], cons(s), ics)
+                except NotImplementedError:
+                    continue
                 rv1.append(s.subs(solved_constants))
+            if len(rv1) == 1:
+                return rv1[0]
             rv = rv1
     return rv
 
@@ -796,7 +801,10 @@ def solve_ics(sols, funcs, constants, ics):
                 sol2 = sol
                 sol2 = sol2.subs(x, x0)
                 sol2 = sol2.subs(funcarg, value)
-                subs_sols.append(sol2)
+                # This check is necessary because of issue #15724
+                if not isinstance(sol2, BooleanAtom) or not subs_sols:
+                    subs_sols = [s for s in subs_sols if not isinstance(s, BooleanAtom)]
+                    subs_sols.append(sol2)
 
     # TODO: Use solveset here
     try:
@@ -817,9 +825,6 @@ def solve_ics(sols, funcs, constants, ics):
     if len(solved_constants) > 1:
         raise NotImplementedError("Initial conditions produced too many solutions for constants")
 
-    if len(solved_constants[0]) != len(constants):
-        raise ValueError("Initial conditions did not produce a solution for all constants. Perhaps they are under-specified.")
-
     return solved_constants[0]
 
 def classify_ode(eq, func=None, dict=False, ics=None, **kwargs):
diff --git a/sympy/solvers/tests/test_ode.py b/sympy/solvers/tests/test_ode.py
index 9427a7275d6c..8efaaedf2a54 100644
--- a/sympy/solvers/tests/test_ode.py
+++ b/sympy/solvers/tests/test_ode.py
@@ -1075,6 +1075,8 @@ def test_classify_sysode():
 
 def test_solve_ics():
     # Basic tests that things work from dsolve.
+    assert dsolve(f(x).diff(x) - 1/f(x), f(x), ics={f(1): 2}) == \
+        Eq(f(x), sqrt(2 * x + 2))
     assert dsolve(f(x).diff(x) - f(x), f(x), ics={f(0): 1}) == Eq(f(x), exp(x))
     assert dsolve(f(x).diff(x) - f(x), f(x), ics={f(x).diff(x).subs(x, 0): 1}) == Eq(f(x), exp(x))
     assert dsolve(f(x).diff(x, x) + f(x), f(x), ics={f(0): 1,
@@ -1103,11 +1105,29 @@ def test_solve_ics():
     assert solve_ics([Eq(f(x), C1*sin(x) + C2*cos(x))], [f(x)], [C1, C2],
         {f(0): 1, f(x).diff(x).subs(x, 0): 1}) == {C1: 1, C2: 1}
 
-    # XXX: Ought to be ValueError
-    raises(NotImplementedError, lambda: solve_ics([Eq(f(x), C1*sin(x) + C2*cos(x))], [f(x)], [C1, C2], {f(0): 1, f(pi): 1}))
+    assert solve_ics([Eq(f(x), C1*sin(x) + C2*cos(x))], [f(x)], [C1, C2], {f(0): 1}) == \
+        {C2: 1}
+
+    # Some more complicated tests Refer to PR #16098
+
+    assert dsolve(f(x).diff(x)*(f(x).diff(x, 2)-x), ics={f(0):0, f(x).diff(x).subs(x, 1):0}) == \
+        [Eq(f(x), 0), Eq(f(x), x ** 3 / 6 - x / 2)]
+    assert dsolve(f(x).diff(x)*(f(x).diff(x, 2)-x), ics={f(0):0}) == \
+        [Eq(f(x), 0), Eq(f(x), C2*x + x**3/6)]
+
+    K, r, f0 = symbols('K r f0')
+    sol = Eq(f(x), -K*f0*exp(r*x)/((K - f0)*(-f0*exp(r*x)/(K - f0) - 1)))
+    assert (dsolve(Eq(f(x).diff(x), r * f(x) * (1 - f(x) / K)), f(x), ics={f(0): f0})) == sol
+
+
+    #Order dependent issues Refer to PR #16098
+    assert dsolve(f(x).diff(x)*(f(x).diff(x, 2)-x), ics={f(x).diff(x).subs(x,0):0, f(0):0}) == \
+        [Eq(f(x), 0), Eq(f(x), x ** 3 / 6)]
+    assert dsolve(f(x).diff(x)*(f(x).diff(x, 2)-x), ics={f(0):0, f(x).diff(x).subs(x,0):0}) == \
+        [Eq(f(x), 0), Eq(f(x), x ** 3 / 6)]
 
     # XXX: Ought to be ValueError
-    raises(ValueError, lambda: solve_ics([Eq(f(x), C1*sin(x) + C2*cos(x))], [f(x)], [C1, C2], {f(0): 1}))
+    raises(NotImplementedError, lambda: solve_ics([Eq(f(x), C1*sin(x) + C2*cos(x))], [f(x)], [C1, C2], {f(0): 1, f(pi): 1}))
 
     # Degenerate case. f'(0) is identically 0.
     raises(ValueError, lambda: solve_ics([Eq(f(x), sqrt(C1 - x**2))], [f(x)], [C1], {f(x).diff(x).subs(x, 0): 0}))


Code review comments:
diff_hunk: @@ -1103,11 +1105,29 @@ def test_solve_ics():
     assert solve_ics([Eq(f(x), C1*sin(x) + C2*cos(x))], [f(x)], [C1, C2],
         {f(0): 1, f(x).diff(x).subs(x, 0): 1}) == {C1: 1, C2: 1}
 
-    # XXX: Ought to be ValueError
-    raises(NotImplementedError, lambda: solve_ics([Eq(f(x), C1*sin(x) + C2*cos(x))], [f(x)], [C1, C2], {f(0): 1, f(pi): 1}))
+    assert solve_ics([Eq(f(x), C1*sin(x) + C2*cos(x))], [f(x)], [C1, C2], {f(0): 1}) == \
+        {C2: 1}
+
+    # Some more complicated tests Refer to PR #16098
+
+    assert dsolve(f(x).diff(x)*(f(x).diff(x, 2)-x), ics={f(0):0, f(x).diff(x).subs(x, 1):0}) == \
+        [Eq(f(x), 0), Eq(f(x), x ** 3 / 6 - x / 2)]
+    assert dsolve(f(x).diff(x)*(f(x).diff(x, 2)-x), ics={f(0):0}) == \
+        [Eq(f(x), 0), Eq(f(x), C2*x + x**3/6)]
+
+    K, r, f0 = symbols('K r f0')
+    sol = Eq(f(x), -K*f0*exp(r*x)/((K - f0)*(-f0*exp(r*x)/(K - f0) - 1)))
+    assert (dsolve(Eq(f(x).diff(x), r * f(x) * (1 - f(x) / K)), f(x), ics={f(0): f0})) == sol
+
+
+    #Order dependent issues Refer to PR #16098
+    assert dsolve(f(x).diff(x)*(f(x).diff(x, 2)-x), ics={f(x).diff(x).subs(x,0):0, f(0):0}) == \
+        [Eq(f(x), 0), Eq(f(x), x ** 3 / 6)]
+    assert dsolve(f(x).diff(x)*(f(x).diff(x, 2)-x), ics={f(0):0, f(x).diff(x).subs(x,0):0}) == \
+        [Eq(f(x), 0), Eq(f(x), x ** 3 / 6)]
 
     # XXX: Ought to be ValueError
-    raises(ValueError, lambda: solve_ics([Eq(f(x), C1*sin(x) + C2*cos(x))], [f(x)], [C1, C2], {f(0): 1}))
+    raises(NotImplementedError, lambda: solve_ics([Eq(f(x), C1*sin(x) + C2*cos(x))], [f(x)], [C1, C2], {f(0): 1, f(pi): 1}))
path: sympy/solvers/tests/test_ode.py
text: @user1:
Okay can we change it so that this NotImplementedError becomes ValueError?

@author:
It was meant to be `NotImplementedError` at the first place. Should I change it to `ValueError`

@user1:
Yes make it so that this example raises ValueError.