Solve motion equations for first ODE using scipy

Question

I would like to solve motion first order ODE equations using scipy solve_ivp function. I can see that I'm doing something wrong because this should be an ellipse but I'm plotting only four points. Are you able to spot the mistake?

import math
import matplotlib.pyplot as plt 
import numpy as np 
import scipy.integrate


gim = 4*(math.pi**2)
x0 = 1 #x-position of the center or h
y0 = 0 #y-position of the center or k
vx0 = 0 #vx position
vy0 = 1.1* 2* math.pi #vy position
initial = [x0, y0, vx0, vy0] #initial state of the system
time = np.arange(0, 1000, 0.01) #period

def motion(t, Z): 
  dx = Z[2] # vx
  dy = Z[3] # vy
  dvx = -gim/(x**2+y**2)**(3/2) * x * Z[2]
  dvy = -gim/(x**2+y**2)**(3/2) * y * Z[3]
  return [dx, dy, dvx, dvy]

sol = scipy.integrate.solve_ivp(motion, t_span=time, y0= initial, method='RK45')
plt.plot(sol.y[0],sol.y[1],"x", label="Scipy RK45 solution")
plt.show()

I don't think you've posted the exact code you're using: `x` and `y` are not defined in `motion`. `solve_ivp` argument `t_span` expects the end points of the time interval (i.e., two numbers). Finally, what is this modelling? Having dvx depend on vx (and similarly for dvy) introduces damping, which I don't *think* will produce an ellipse; the equations are nonlinear though, so it's difficult to say. — Rory Yorke, Nov 20 '21 at 15:42

score 0 · Answer 1 · answered Nov 20 '21 at 16:02

The code should not be able to run from a fresh workspace. The variables x,y that are used in the gravitation formula are not declared anywhere. So insert the line x,y = Z[:2] or similar.

The gravitation formula usually does not contain the velocity components. Remove Z[2], Z[3].

Check again what the time span and evaluation times arguments expect. The time span takes the first two values from the array. So change to t_span=time[[0,-1]] to build the pair of first and last time value.

The second plot suffers from insufficient evaluation points, the line segments used are too large. With your time array that should not be a problem.

Solve motion equations for first ODE using scipy

1 Answers1