dsolve fails: Initial conditions produced too many solutions for constants

[oscarbenjamin]

For this ODE dsolve can find the general solution but fails when asked to solve with initial conditions:

In [1]: eq = diff(f(x), x) - sqrt(f(x)**2 - 1)

In [2]: eq
Out[2]: 
     ___________           
    ╱  2           d       
- ╲╱  f (x) - 1  + ──(f(x))
                   dx      

In [3]: dsolve(eq)
Out[3]: f(x) = cosh(C₁ + x)

In [5]: dsolve(eq, ics={f(0):2})
...
NotImplementedError: Initial conditions produced too many solutions for constants

What happens here is that it tried to solve for C1 like:

In [6]: gensol = dsolve(eq)

In [7]: gensol
Out[7]: f(x) = cosh(C₁ + x)

In [8]: gensol.subs(x, 0).subs(f(0), 2)
Out[8]: 2 = cosh(C₁)

In [9]: solve(gensol.subs(x, 0).subs(f(0), 2), Symbol('C1'))
Out[9]: [log(2 - √3), log(√3 + 2)]

We find that the expression for the general solution can satisfy the initial condition for two different values of the constant $C_1$. This is because there are two different intersecting cosh curves crossing through the f(0)=2 point in the x,f plane like:

While both curves satisfy the initial condition they have different slopes there and only one of them satisfies the ODE at that point which requires that the gradient be +sqrt(3) there:

In [14]: eq
Out[14]: 
     ___________           
    ╱  2           d       
- ╲╱  f (x) - 1  + ──(f(x))
                   dx      

In [15]: solve(eq, f(x).diff(x))
Out[15]: 
⎡   ___________⎤
⎢  ╱  2        ⎥
⎣╲╱  f (x) - 1 ⎦

In [16]: solve(eq, f(x).diff(x))[0].subs(x, 0).subs(f(0), 2)
Out[16]: √3

In general it would be better if dsolve passed both the ODE and the initial conditions when solving for the constants:

In [20]: sys = [eq.subs(f(x), gensol.rhs).doit().subs(x, 0), gensol.subs(x, 0).subs(f(0), 2)]

In [21]: sys
Out[21]: 
⎡     _______________                         ⎤
⎢    ╱     2                                  ⎥
⎣- ╲╱  cosh (C₁) - 1  + sinh(C₁), 2 = cosh(C₁)⎦

In [22]: solve(sys, Symbol('C1'))
Out[22]: [(log(√3 + 2),)]

This gives the correct values for $C_1$ that satifies both the ODE and the initial condition at x=0.

In this particular case the ODE is solved as a separable equation which initially obtains an implicit equation for f(x):

In [23]: classify_ode(eq)
Out[23]: 
('separable',
 '1st_exact',
 '1st_power_series',
 'lie_group',
 'separable_Integral',
 '1st_exact_Integral')

In [24]: dsolve(eq, hint='separable_Integral')
Out[24]: 
f(x)                             
 ⌠                               
 ⎮        1                ⌠     
 ⎮   ─────────── dy = C₁ + ⎮ 1 dx
 ⎮      ________           ⌡     
 ⎮     ╱  2                      
 ⎮   ╲╱  y  - 1                  
 ⌡                               
                                 

In [25]: dsolve(eq, hint='separable_Integral').doit()
Out[25]: 
   ⎛   ___________       ⎞         
   ⎜  ╱  2               ⎟         
log⎝╲╱  f (x) - 1  + f(x)⎠ = C₁ + x

In a pen and paper calculation you would definitely apply the initial condition to this implicit equation rather than solving for f(x) first and then trying to find C1 later:

In [26]: eq_ics = dsolve(eq, hint='separable_Integral').doit().subs(x, 0).subs(f(0), 2)

In [27]: eq_ics
Out[27]: log(√3 + 2) = C₁

In [28]: eq_implicit = dsolve(eq, hint='separable_Integral').doit().subs(eq_ics.rhs, eq_ics.lhs)

In [29]: eq_implicit
Out[29]: 
   ⎛   ___________       ⎞                  
   ⎜  ╱  2               ⎟                  
log⎝╲╱  f (x) - 1  + f(x)⎠ = x + log(√3 + 2)

In [30]: solve(eq_implicit, f(x))
Out[30]: [cosh(x + log(√3 + 2))]

There are two things that can be improved here:

1. Apply initial conditions to the implicit equation.
2. Also require solved values of the constants to satisfy the ODE.

Another reason why dsolve might fail like this is if the initial value problem genuinely does not have a unique solution. Currently nothing is done to check for that case but potentially it would make sense for dsolve to return multiple solutions then.