odeint(deriv_x, xinit, t)
uses xinit
as its initial guess for x
. This value for x
is used when evaluating deriv_x
.
deriv_x(xinit, t)
raises a divide-by-zero error since x[0] = xinit[0]
equals 0, and deriv_x
divides by x[0]
.
It looks like you are trying to solve the second-order ODE
r'' = - C rhat
---------
|r|**2
where rhat is the unit vector in the radial direction.
You appear to be separating the x
and y
coordinates into separate second-order ODES:
x'' = - C y'' = - C
----- and -----
x**2 y**2
with initial conditions x0 = 0 and y0 = 20056.
This is very problematic. Among the problems is that when x0 = 0
, x''
blows up. The original second-order ODE for r''
does not have this problem -- the denominator does not blow up when x0 = 0
because y0 = 20056
, and so r0 = (x**2+y**2)**(1/2)
is far from zero.
Conclusion: Your method of separating the r''
ODE into two ODEs for x''
and y''
is incorrect.
Try searching for a different way to solve the r''
ODE.
Hint:
- What if your "state" vector is
z = [x, y, x', y']
? - Can you write down a first-order ODE for
z'
in terms ofx
,y
,x'
andy'
? - Can you solve it with one call to
integrate.odeint
?