Ask Your Question

coupled second order differential equation

asked 2013-02-06 10:19:16 +0200

anonymous user


updated 2013-02-06 11:46:20 +0200

calc314 gravatar image


I am working on the restricted three body problem, and it turns out it is a coupled second ODE.


a,t = var('a,t')
x = function('x', t)
y = function('y', t)
dx2 = diff(x,t,2) == -(1-a)*((x-a)/(sqrt((x-a)^2 + y^2))^3)-a*((x+1-a)/(sqrt((x+1-a)^2 + y^2))^3)+x+2*diff(y,t,1)
dy2 = diff(y,t,2) == -(1-a)*(y/(sqrt((x-a)^2 + y^2))^3)-a*(y/(sqrt((x+1-a)^2 + y^2))^3)+y-2*diff(x,t,1)
desolve(dx2, x, ics=[0,1])


ValueError Traceback (most recent call last)

/home/eric/<ipython console> in <module>()

/home/eric/programs/sage-5.5/local/lib/python2.7/site-packages/sage/calculus/desolvers.pyc in desolve(de, dvar, ics, ivar, show_method, contrib_ode)
418 ivars = [t for t in ivars if t is not dvar]
419 if len(ivars) != 1:
--> 420 raise ValueError, "Unable to determine independent variable, please specify."
421 ivar = ivars[0]
422 def sanitize_var(exprs):

ValueError: Unable to determine independent variable, please specify.

I am not sure whats going on here. From what I have looked up on how to solve ODE in sage my code should be correct. I am using cantor for a graphical interface. I have also tried the code in the command line, and still get the same error.

edit retag flag offensive close merge delete

2 Answers

Sort by ยป oldest newest most voted

answered 2013-02-06 15:43:52 +0200

calc314 gravatar image

updated 2013-02-06 21:31:12 +0200

If you want a numerical solution, you can try the following.

dx1 = u
du = -(1-a)*((x-a)/(sqrt((x-a)^2 + y^2))^3)-a*((x+1-a)/(sqrt((x+1-a)^2 + y^2))^3)+x+2*v
dy1 = v
dv = -(1-a)*(y/(sqrt((x-a)^2 + y^2))^3)-a*(y/(sqrt((x+1-a)^2 + y^2))^3)+y-2*u

html('t vs. x,y')
html('phase plane (x vs. y)')
html('t vs. u,v')

fu(x,u,y,v) = -(1-a)*((x-a)/(sqrt((x-a)^2 + y^2))^3)-a*((x+1-a)/(sqrt((x+1-a)^2 + y^2))^3)+x+2*v
fv(x,u,y,v) = -(1-a)*(y/(sqrt((x-a)^2 + y^2))^3)-a*(y/(sqrt((x+1-a)^2 + y^2))^3)+y-2*u

tudot=[(tspace[i],fu(result[i][0],result[i][1],result[i][2],result[i][3])) for i in range(0,len(tspace))]

tvdot=[(tspace[i],fv(result[i][0],result[i][1],result[i][2],result[i][3])) for i in range(0,len(tspace))]


You can run the code at the single cell server by clicking here.

edit flag offensive delete link more


thank you that will work as well.

eric.c.kangas gravatar imageeric.c.kangas ( 2013-02-06 16:26:56 +0200 )edit

Also I need to figure out how to solve numerically the acceleration. I tried also solving for udot, and vdot, but I get an error. error: ValueError Traceback (most recent call last) /home/eric/programs/sage-5.5/local/lib/python2.7/site-packages/sage/calculus/desolvers.pyc in func(y, t) 1538 v = list(y[:]) 1539 v.append(t) -> 1540 return [dec(*v) for dec in desc] 1541 1542 if not compute_jac: /home/eric/programs/sage-5.5/local/lib/python2.7/site-packages/sage/ext/interpreters/ in sage.ext.interpreters.wrapper_rdf.Wrapper_rdf.__call__ (sage/ext/interpreters/wrapper_rdf.c:1482)() 73 74 ---> 75 76 77 ValueError: --------------------------------------------------------------------------- error Traceback (most recent call last)

eric.c.kangas gravatar imageeric.c.kangas ( 2013-02-06 18:39:56 +0200 )edit

I've edited the answer above to include a straight-forward computation of the acceleration from the velocity curves.

calc314 gravatar imagecalc314 ( 2013-02-06 19:08:05 +0200 )edit

interesting so the odeint function isn't setup yet to calculate acceleration at this time.

eric.c.kangas gravatar imageeric.c.kangas ( 2013-02-06 19:40:39 +0200 )edit

Actually, we don't need odeint to compute the acceleration, and I didn't need to do the difference quotient either. We just need to use the ODEs for the 2nd derivative, which is the acceleration. See my latest edit above.

calc314 gravatar imagecalc314 ( 2013-02-06 21:29:18 +0200 )edit

answered 2013-02-06 12:05:49 +0200

calc314 gravatar image

updated 2013-02-06 12:13:47 +0200

You'll need to let Sage know what the independent variable is by including ivar=t in your command. Further, since you are solving a system of equations, you need to use desolve_system.

a,t = var('a,t')
x = function('x', t)
y = function('y', t)
dx2 = 0 == -diff(x,t,2)-(1-a)*((x-a)/(sqrt((x-a)^2 + y^2))^3)-a*((x+1-a)/(sqrt((x+1-a)^2 + y^2))^3)+x+2*diff(y,t,1)
dy2 = 0 == -diff(y,t,2)-(1-a)*(y/(sqrt((x-a)^2 + y^2))^3)-a*(y/(sqrt((x+1-a)^2 + y^2))^3)+y-2*diff(x,t,1)

This code gets farther but gives an error as well. This might be related to limitations in Maxima's inverse Laplace transform abilities.

edit flag offensive delete link more


I will try out scipy

eric.c.kangas gravatar imageeric.c.kangas ( 2013-02-06 13:31:21 +0200 )edit

That might work well. Note that you can run scipy from within Sage.

calc314 gravatar imagecalc314 ( 2013-02-06 14:44:47 +0200 )edit

okay even I am having issues with scipy. I try importing odeint but it won't let me. code: import scipy; from scipy import odeint; error: --------------------------------------------------------------------------- ImportError Traceback (most recent call last) /home/eric/<ipython console=""> in <module>() ImportError: cannot import name odeint

eric.c.kangas gravatar imageeric.c.kangas ( 2013-02-06 16:25:19 +0200 )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


Asked: 2013-02-06 10:19:16 +0200

Seen: 2,625 times

Last updated: Feb 06 '13