sympy / sympy

A computer algebra system written in pure Python
https://sympy.org/
Other
12.79k stars 4.39k forks source link

SymPy gives a wrong solution for an ODE of first order #26593

Closed terben closed 4 months ago

terben commented 4 months ago

I am trying to solve with SymPy the ODEs:

$$- \sqrt{y_1^{2}{\left(x \right)} - 1} + \frac{d}{d x} y_1{\left(x \right)} = 0$$

and

$$- \sqrt{\left(y_2{\left(x \right)} + 1\right)^{2} - 1} + \frac{d}{d x} y_2{\left(x \right)}=0$$

The solution for the first ODE is $y_1(x)=\cosh(C+x)$ which I can reproduce with SymPy. The solution for the second ODE is $y_2(x)=\cosh(C+x)-1$. This can be seen with the substitution $y_3(x)=y_2(x)+1$. It produces for $y_3(x)$ the same ODE as for $y_1(x)$.

However, SymPy yields the following solution for $y_2(x)$:

$$y2{\left(x \right)} = e^{- C{1} - x} + \frac{e^{C_{1} + x}}{4} - 1$$

which cannot be transformed into $y_2(x)=\cosh(C_1+x)-1$.

Here is my interaction with SymPy:

In [5]: x = sp.symbols('x', real=True)

In [6]: y = sp.symbols('y', cls=sp.Function)

In [7]: y = y(x)

In [8]: dgl1 = sp.diff(y, x) - sp.sqrt((y)**2 - 1)

In [9]: sp.dsolve(dgl1)
Out[9]: Eq(y(x), cosh(C1 + x))

In [10]: dgl2 = sp.diff(y, x) - sp.sqrt((y + 1)**2 - 1)

In [11]: sp.dsolve(dgl2)
Out[11]: Eq(y(x), exp(-C1 - x) + exp(C1 + x)/4 - 1)

In [14]: sp.__version__
Out[14]: '1.12'
oscarbenjamin commented 4 months ago

The solution is equivalent:

In [438]: sol = dsolve(dgl2)

In [439]: sol
Out[439]:
                   C₁ + x
        -C₁ - x   ℯ
y(x) = ℯ        + ─────── - 1
                     4

In [440]: C1, C2 = symbols('C1, C2')

In [441]: sol.subs(C1, C2+log(2))
Out[441]:
        -C₂ - x    C₂ + x
       ℯ          ℯ
y(x) = ──────── + ─────── - 1
          2          2

In [442]: sol.subs(C1, C2+log(2)).rewrite(cos).simplify()
Out[442]: y(x) = cosh(C₂ + x) - 1
terben commented 4 months ago

Oh yes. Thanks a lot! I did not see this despite trying for quite some time :-)

Of course it would be nice if SymPy gave immediately the cosh-expression but of course this is only an inconvenience and not an error. From my side we can close this issue.

asmeurer commented 4 months ago

That's an issue for the integrator, not dsolve. Although even that is not a clear cut thing because some people may prefer for the integrator to return hyperbolic trig functions when possible and some may prefer to always get exponentials.

oscarbenjamin commented 4 months ago

I think that simplify is being used somewhere and the difference is:

In [46]: simplify(exp(x) + exp(-x))
Out[46]: 2⋅cosh(x)

In [47]: simplify(exp(x) + exp(-x)/4)
Out[47]: 
      -x
 x   ℯ  
ℯ  + ───
      4 

The integrator gives an answer with log:

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

In [54]: eq2 = dsolve(eq, hint='separable_Integral').doit()

In [55]: eq2
Out[55]: 
   ⎛   ___________       ⎞         
   ⎜  ╱  2               ⎟         
log⎝╲╱  f (x) - 1  + f(x)⎠ = C₁ + x

In [56]: solve(eq2, f(x))
Out[56]: [cosh(C₁ + x)]

In [57]: solve(eq2, f(x), simplify=False)
Out[57]: 
⎡ -C₁ - x    C₁ + x⎤
⎢ℯ          ℯ      ⎥
⎢──────── + ───────⎥
⎣   2          2   ⎦
oscarbenjamin commented 4 months ago

I think we can close this.

I've opened gh-26595 as a followup for a different problem with this ODE.