# desolve_system_rk4 vs desolve_odeint

I am getting different results (a result vs no result) from desolve_system_rk4 and desolve_odeint for the same system. What am I missing? My code is the following:

reset()
from sage.calculus.desolvers import desolve_system_rk4,desolve_odeint
var('eta,r,ph');

drdt=-sqrt(2.0*r + 200.0 / r - 100.0)*(r - 2.0) / (r^2)
dphdt=10.0*(r - 2.0) / (r^3)

sol1=desolve_system_rk4([drdt,dphdt],[r,ph],ics=[0,4.0,0.3],ivar=eta,end_points=1000,step=0.01)
sol2=desolve_odeint([drdt,dphdt],[4.0,0.3],srange(0,1000,0.01),[r,ph])

D=[[j*cos(k),j*sin(k)] for i,j,k in sol1]

p1=list_plot(D,aspect_ratio=1)
p2=line(zip(sol2[:,0]*cos(sol2[:,1]),sol2[:,0]*sin(sol2[:,1])))

show(graphics_array([p1,p2]))

edit retag close merge delete

Sort by ยป oldest newest most voted

For your initial conditions, drdt has an imaginary part. Differential equation solvers in SageMath generally only work with real values, so desolve_odeint using SciPy doesn't give an answer as expected. The solver desolve_system_rk4 uses Maxima and shouldn't really be returning an answer: that looks like a Maxima bug.

If you set r(0)=1 to keep the system real, you get the same answer from both solvers.

more

Thank you Paul. It is a more serious issue than I thought. I will try to reformulate my physical problem.

( 2017-01-07 07:52:19 +0200 )edit