2
$\begingroup$

I have written a program to solve for Newton's 2nd Law of motion for a given force law, in 2D polar coordinates. It is known that if the force law is of the form $k/r^2$,we get conic sections as solution(https://en.wikipedia.org/wiki/Kepler_problem).

But when I run the program, I get a spiral trajectory (which isn't one of the conic sections).

Componentwise, Newton's 2nd law for central force $F(r)=k/r^2$(and unit mass) are-$$\frac{d^2r}{dt^2}=k/r^2$$ $$\frac{d^2\theta}{dt^2}=0$$ where the first equation is the radial part and the 2nd equation is for angular part(which vanishes for Central force problems,due to conservation of angular momentum).

The equations to be used in the program are a set of 4 coupled first order ODEs, derived from the above 2nd order ODEs.(Since odeint cant directly handle 2nd order ODE).The equations I used in the program were--- $$dr/dt=v_r$$$$dv_r/dt=k/r^2$$$$d\theta/dt=v_\theta$$$$dv_\theta/dt=0$$with $k=1$ for simplicity.

I have tried with many possible initial conditions, and tried running for different values of $k$ to ensure. Why is this happening? Is this because of the choice of initial condition or some faulty implementation or some numerical problem? Or is their a problem with the physics?

Below is my minimal code

def vec(w, t):
    r, vr, theta, vtheta = w
    return [vr, r**(-2.0), vtheta, 0]

def newton(vec, initial, t):
    wsol = odeint(vec, initial, t)
    return [wsol[:, 0], wsol[:, 2]]

T = np.linspace(0, 50, 1000)    
initial = [2, 10, radians(0),2] #The order is radius, radial velocity, angle, angular velocity
R = newton(vec, initial, T)[0]
Theta = newton(vec, initial, T)[1]
plt.polar(Theta, R, "r", lw="1")
plt.show()
$\endgroup$
1
  • $\begingroup$ I couldn't find much good tags suitable for this question. People who are more familiar with the tags, please help me edit the question with more suitable tags. $\endgroup$ Jul 27, 2020 at 8:47

1 Answer 1

4
$\begingroup$

The correct dynamic equations in the polar coordinates should be

$ \dot{v_r} = \omega^2 r - \alpha/r^2 \\ \dot{\omega} = - 2 v_r \omega /r\\ \dot{\theta} = \omega \\ \dot{r} = v_r $

Here is the fixed Python code:

from math import *
import numpy as np
from matplotlib import pyplot as plt
from scipy.integrate import odeint


def vec(w,t):
    r,vr,theta,omega=w
    return [vr, r*omega**2 - 1e0/r**2, omega, -2*vr*omega/r]

def newton(vec,initial,t):
    wsol=odeint(vec,initial,t)
    return [wsol[:,0],wsol[:,1],wsol[:,2],wsol[:,3]]

t=np.linspace(0,4,1121)
#The order is radius,radial velocity,angle,angular velocity.
initial=[0.5, -1e0, radians(0), 1e0] 

r=newton(vec,initial,t)[0]
vr=newton(vec,initial,t)[1]
theta=newton(vec,initial,t)[2]
omega=newton(vec,initial,t)[3]

plt.polar(theta,r, "r", lw="1")
plt.show()

Here is the output:

enter image description here

$\endgroup$
1
  • $\begingroup$ Oh! I used zero angular velocity by mistake. I forgot that angular acceleration isn't 0 when distance from axis of rotation changes.Got it.Thank you. $\endgroup$ Jul 27, 2020 at 16:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.