Revision history [back]

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

followed by:

line(tx)+line(ty,color='green')

and

line(xy)

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

followed by:

line(tx)+line(ty,color='green')

and

line(xy)

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


udot=[]
vdot=[]
u=result[:,1]
v=result[:,3]
for i in range(0,len(tspace)-1):
    udot.append((u[i+1]-u[i])/(tspace[i+1]-tspace[i]))
    vdot.append((v[i+1]-v[i])/(tspace[i+1]-tspace[i]))
tudot=zip(tspace,udot)
tvdot=zip(tspace,vdot)

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

line(tx)+line(ty,color='green')

line(xy)

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

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


udot=[]
vdot=[]
u=result[:,1]
v=result[:,3]
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)-1):
    udot.append((u[i+1]-u[i])/(tspace[i+1]-tspace[i]))
    vdot.append((v[i+1]-v[i])/(tspace[i+1]-tspace[i]))
tudot=zip(tspace,udot)
tvdot=zip(tspace,vdot)
range(0,len(tspace))]
 html('t vs. accelerations')
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.