# SymPy to Solve Initial Value Condition (First Order Differential Equation)

Hi all,

I want to check is my code correct to compute IVP of

\frac{dx}{dt} = \sqrt{2g} \frac{R}{\sqrt{R+x}}

the initial condition is x(0)=0

using SymPy
@syms t, g, R
@syms x()

# Computing IVP with SymPy
x(t).diff(t)
diffeq = Eq(x(t).diff(t), sqrt(2g)*R*(1/sqrt(R+x(t))) ); string(diffeq)

# To solve the ODE, pass it and the function to solve for to dsolve.
xt = dsolve(diffeq, x(t))
diffeq3 = xt[3](Dict(R => 4000, g => 78545 ))
# To solve the Initial Value Problem, initial condition is 0
ivp = dsolve(diffeq3, x(t), ics = Dict(x(0) => 0 ))


there are 3 solutions at xt and when I try to compute IVP it becomes error:

PyError (\$(Expr(:escape, :(ccall(#= /home/browni/.julia/packages/PyCall/ygXW2/src/pyfncall.jl:43 =# @pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, pyargsptr, kw))))) <class âValueErrorâ>
ValueError(âCouldnât solve for initial conditionsâ)
** File â/home/browni/.julia/conda/3/envs/lasthrim_env/lib/python3.10/site-packages/sympy/solvers/ode/ode.pyâ, line 640, in dsolve**
** return _helper_simplify(eq, hint, hints, simplify, ics=ics)**
** File â/home/browni/.julia/conda/3/envs/lasthrim_env/lib/python3.10/site-packages/sympy/solvers/ode/ode.pyâ, line 695, in _helper_simplify**
** solved_constants = solve_ics([rv], [r[âfuncâ]], cons(rv), ics)**
** File â/home/browni/.julia/conda/3/envs/lasthrim_env/lib/python3.10/site-packages/sympy/solvers/ode/ode.pyâ, line 802, in solve_ics**
** raise ValueError(âCouldnât solve for initial conditionsâ)**

Stacktrace:

This isnât quite the way to set this up. You are calling dsolve twice, not once. However, here SymPy also gives up. You can do the work manually, though it isnât so pretty (sorry, I switched some variable names):

julia> eq = diff(u(t), t) - â(2g) * R / â(R + u(t))

julia> sols =  dsolve(eq);

julia> sol = rhs(sols[3]);

julia> c1 = free_symbols(sol)[1]
Câ

julia> C1 = solve(sol(t=>0) ~ 0, c1)[2];

julia> sol(c1 => C1)
...

1 Like

Yes thanks @j_verzani , I am starting to learn how to use SymPy to solve DE now.

this is working code:

using SymPy
@syms t, g, R
@syms x()

# Computing IVP with SymPy
diffeq = diff(x(t), t) - â(2g) * R / â(R + x(t))

# To solve the ODE, pass it and the function to solve for to dsolve.
xt = dsolve(diffeq)

# To solve the Initial Value Problem, initial condition is 0
sol_xt = rhs(xt[3])
c1 = free_symbols(sol_xt)[1]
t1 = free_symbols(sol_xt)[4]

C1 = solve(sol_xt(t=>0) ~ 0, c1)[2]
T1 = solve(sol_xt(c1=>C1) + 4000 ~ 240000, t1)[2]

ivp = sol_xt(c1=>C1)
#ivp = sol_xt(c1=>C1,R => 4000, g => 78545)

println("Modeling with First Order Equations with Julia")
println("A Rocket is Launched straight up from the surface of the earth")
println("Neglect air resistance")
println("The initial condition: x(0) = 0")
println("g = GM/(R^2)")
println("g = 78,545 mi/hr^2")

println("")
println("The time required for the rocket to go 240,000 miles (in hours): ")
println(T1(R => 4000, g => 78545).evalf())

println("")
println("The initial value problem solution is:")
println("x(t) = ")
ivp