# coupled second order differential equation Anonymous

Hi,

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

code:

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


error:

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

Sort by » oldest newest most voted

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

a=0.5
var('x,y,u,v,t')
tspace=srange(0,10,.05)
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
result=desolve_odeint([dx1,du,dy1,dv],[1,0,1,0],tspace,[x,u,y,v],ivar=t)
xy=zip(result[:,0],result[:,2])
tx=zip(tspace,result[:,0])
ty=zip(tspace,result[:,2])
tu=zip(tspace,result[:,1])
tv=zip(tspace,result[:,3])

html('t vs. x,y')
show(line(tx)+line(ty,color='green'),axes_labels=['$t$','$x,y$'])
html('phase plane (x vs. y)')
show(line(xy),axes_labels=['$x$','$y$'],aspect_ratio=1)
html('t vs. u,v')
show(line(tu)+line(tv,color='green'),axes_labels=['$t$','$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],result[i],result[i],result[i])) for i in range(0,len(tspace))]

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

show(line(tudot,color='red')+line(tvdot,color='orange'),axes_labels=['$t$','$\dot{u},\dot{v}$'])


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

more

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/wrapper_rdf.so 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)

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

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

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.

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)
desolve_system([dx2,dy2],[x,y],ivar=t)


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

more

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

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