Ask Your Question
0

desolve_system_rk4 vs desolve_odeint

asked 2017-01-06 14:36:03 +0100

anonymous user

Anonymous

updated 2017-01-06 20:23:29 +0100

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 flag offensive close merge delete

1 Answer

Sort by ยป oldest newest most voted
0

answered 2017-01-06 22:18:33 +0100

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.

edit flag offensive delete link more

Comments

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

implicitnone gravatar imageimplicitnone ( 2017-01-07 07:52:19 +0100 )edit

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: 2017-01-06 14:29:26 +0100

Seen: 434 times

Last updated: Jan 06 '17