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

General Usage
1

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:

2

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
3

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