3

I have the following differential equation that I would like to solve with SymPy

enter image description here

This differential equation has the implicit solution (with h(0) = [0,1) and t = [0, inf) )

enter image description here

but SymPy gives

enter image description here

which other packages such as Maxima are able to find. With SymPy I am unable, however. Is there a way to do so? My code is

import sympy as sp
sp.init_printing(use_unicode=True)
h = sp.symbols('h', function=True)
t = sp.symbols('t')
eq = sp.Eq(sp.Derivative(h(t),t), (1 - h(t))**sp.Rational(4,3) / h(t))
sp.dsolve(eq)
Chiel
  • 6,006
  • 2
  • 32
  • 57

1 Answers1

4

SymPy leaves the integral unevaluated because it is unsure about the sign of 1-y in the integral.

The differential equation has a singularity at h=1, and its behavior depends on what side of 1 we are. There isn't a way to say that h(t) < 1, but one can substitute h(t) = 1 - g(t) where g is a positive function: 

g = sp.symbols('g', function=True, positive=True)
eq1 = eq.subs(h(t), 1 - g(t))
print(sp.dsolve(eq1))

This returns an explicit solution of the ODE (actually three of them, as SymPy solves a cubic equation). The first one of those looks reasonable.

Eq(g(t), (-2*(C1 + t)/(sqrt(-8*(C1 + t)**3 + 729) + 27)**(1/3) - (sqrt(-8*(C1 + t)**3 + 729) + 27)**(1/3))**3/27)
  • Interesting. Is there not a way to make some assumptions on `h(t)`, i.e. it is smaller than 1 and larger or equal to 0? – Chiel Jul 26 '17 at 00:35
  • "Nonnegative": yes, "less than 1": no. See http://docs.sympy.org/latest/modules/core.html#module-sympy.core.assumptions –  Jul 26 '17 at 00:43
  • Hmm, that is a pity. With this level of complexity, solving by hand seems more straightforward. – Chiel Jul 26 '17 at 11:03