desolve_system_rk4 vs desolve_odeint

Anonymous

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]))

add a comment

1 Answer

Sort by » oldest newest most voted

answered 8 years ago

paulmasson
1580 ●1 ●24 ●45 http://analyticphysics...

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.

link

Comments

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

implicitnone ( 8 years ago )

add a comment

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 8 years ago

Seen: 458 times

Last updated: Jan 06 '17

desolve_system_rk4 vs desolve_odeint savecancel