I have a physics problem and have never used odeint or any numerical solving for ODE on python and am a little confused. I tried looking at other examples but could not understand and am hoping for a little help. My ODE is:
Where α is a given angle and g is a constant.
r=sqrt(x^2+y^2)
The program will ask for x_0, dx/dt_0 and dy/dt_0.
I'm mostly unsure how to solve ODE's in python. I have seen that I should split my ODE into dr'/dt because odeint will only do a first order ODE's. Could someone help explain how to do this?
I tried using another example to do as much as possible but am stuck:
import numpy as np
import matplotlib.pylab as plt
from scipy.integrate import odeint
pi=np.pi
sin=np.sin
cos=np.cos
sqrt=np.sqrt
alpha=pi/4
g=9.80665
y0=0.0
theta0=0.0
x=[]
y=[]
sina = sin(alpha)**2
second_term = g*sin(alpha)*cos(alpha)
x0 = float(raw_input('What is the initial x in meters?'))
x_vel0 = float(raw_input('What is the initial velocity in the x direction in m/s?'))
y_vel0 = float(raw_input('what is the initial velocity in the y direction in m/s?'))
t_f = float(raw_input('What is the maximum time in seconds?'))
r0 = x0
r = sqrt(float(x)**2 + float(y)**2)
def deriv_z(z,r):
r, rdot=z
return [rdot,r*sina-second_term]
zdot0=x_vel0**2+y_vel0**2
z0 = [r0,zdot0]
times = np.linespace(0, t_f, 1000)
z = odeint(deriv_z,z0,times)
There's a great example I found to help me with my planetary orbit ODE solving.
It uses an adaptive step size solver and plots the orbit nicely. If your solution can handle it, you may also try to use the faster 'lsoda' or 'vode' options instead of 'dopri5' (which is a very solid standard).