v=3
t,x,y,dx,dy=var('t,x,y,dx,dy')
r=sqrt((x-v*t)^2+y^2)
f=(tanh(r+4)-tanh(r-4))/(2*tanh(4))
ft=derivative(f,t)
fx=derivative(f,x)
fy=derivative(f,y)
t_tt=f^2*fx*v^3
t_tx=-f*fx*v^2
t_xx=v*fx
t_ty=-(f*fy*v^2)/2
t_xy=fy*v/2
x_tt=f^3*fx*v^4-f*fx*v^2-v*ft
x_tx=-t_tt
x_xx=-t_tx
x_ty=-(f^2*fy*v^3+v*fy)/2
x_xy=(f*fy*v^2)/2
y_tt=-f*fy*v^2
y_tx=v*fy/2
g=-x_tt-x_tx*dx-x_xx*dx^2-x_ty*dy-x_xy*dx*dy+dx*(t_tt+t_tx*dx+t_xx*dx^2+t_ty*dy+t_xy*dx*dy)
h=-y_tt-y_tx*dx+dy*(t_tt+t_tx*dx+t_xx*dx^2+t_ty*dy+t_xy*dx*dy)
z=[dx,dy,g,h]
sol=desolve_odeint(z,[0.1,0.1,0,1],srange(0,10,0.1),[x,y,dx,dy],ivar=t)
p=line(zip(sol[:,0],sol[:,1]))
p.show()
This code gives me an error saying "excess work done," so I'm assuming that the derivatives blew up somewhere, but I don't know why it says this for these initial conditions and not for others. Now, odeint has parameters hmin, hmax, rtol, and atol, but I don't know if there is a way to fix a problem using these.
