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

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

thank you that will work as well.

( 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/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)

( 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.

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

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

( 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.

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

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

I will try out scipy

( 2013-02-06 13:31:21 +0200 )edit

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

( 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

( 2013-02-06 16:25:19 +0200 )edit