Solved ODE by Sympy doesn't satisfy the given initial condition

user_6714669 0
Viewed 93

I'm trying to solve the following simple ODE using dsolve function in Sympy

v = sp.Function('v', real=True)
t, g, m, c = sp.symbols('t g m c', positive=True, real=True)
eq = sp.Eq(v(t).diff(t), g - c/m*v(t)**2)
sol = sp.dsolve(eq, v(t), ics={v(0):0})

The result isn't as simple as expected but the strange point is that if you check the solution to insure that the given initial condition is satisfied, it doesn't return the correct result (v(0)=0), as passed in the dsolve by ics argument.

sol.subs(t,0)
ode sympy
2 Answers

Replacing u=c/m*v and c/m*g=a^2 one gets the equation u'=a^2-u^2 which is separable and leads to

log|u+a|-log|u-a| = 2*a*t+c

This should be transformed to

(u+a)/(u-a) = C*exp(2*a*t)

u = a*(C*exp(2*a*t)+1)/(C*exp(2*a*t)-1)

where then the initial condition leads to C=-1 and u(t)=a*tanh(a*t).

As the constant in the sympy solution is apparently not removed from the exponent, and the resolution of the absolute value towards a negative sign also apparently not followed, sympy ends up with a variation of c=log(-1)/2, which makes no sense in this context of an equation over the reals. Also -1 is on the usual branch cut of the complex logarithm, making this even more uncertain. In the actual solution one finds log((-1)**sqrt(m)) which is not a well-defined expression.

user_30881380

The problem is that dsolve finds the general solution:

In [41]: sol
Out[41]: 
                      √g⋅√m               
v(t) = ───────────────────────────────────
              ⎛             ⎛ C₁⋅√c⋅√g⋅m⎞⎞
              ⎜          log⎝ℯ          ⎠⎟
              ⎜√c⋅√g⋅t - ────────────────⎟
              ⎜                 2        ⎟
       √c⋅tanh⎜──────────────────────────⎟
              ⎝            √m            ⎠

Before attempting to solve for C1 here it would be better to apply some constant simplification:

In [42]: [logterm] = sol.atoms(log)

In [43]: C1 = Symbol('C1')

In [44]: sol2 = sol.subs(logterm, C1)

In [45]: sol2
Out[45]: 
                √g⋅√m         
v(t) = ───────────────────────
              ⎛  C₁          ⎞
              ⎜- ── + √c⋅√g⋅t⎟
              ⎜  2           ⎟
       √c⋅tanh⎜──────────────⎟
              ⎝      √m      ⎠

In [46]: [c] = solve(sol2.rhs.subs(t, 0), C1)

In [47]: c
Out[47]: ⅈ⋅π⋅√m

In [48]: sol2.subs(C1, c)
Out[48]: 
                 √g⋅√m          
v(t) = ─────────────────────────
              ⎛          ⅈ⋅π⋅√m⎞
              ⎜√c⋅√g⋅t - ──────⎟
              ⎜            2   ⎟
       √c⋅tanh⎜────────────────⎟
              ⎝       √m       ⎠

In [49]: sol2.subs(C1, c).simplify()
Out[49]: 
                 ⎛√c⋅√g⋅t⎞
       √g⋅√m⋅tanh⎜───────⎟
                 ⎝   √m  ⎠
v(t) = ───────────────────
                √c

Perhaps actually it is better to solve for the constant at an earlier stage in the solving process. Probably on pen and paper it would make sense to solve for the constant as soon as integration happens and while the equation is still in implicit form.

user_94509910
Related