Why is $1/r^2$ force law giving spiral trajectory?

Computational Science Asked on December 9, 2020

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(

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$$
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$$$$dtheta/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")

One Answer

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):
    return [vr, r*omega**2 - 1e0/r**2, omega, -2*vr*omega/r]

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

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


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

Here is the output:

enter image description here

Correct answer by Maxim Umansky on December 9, 2020

Add your own answers!

Ask a Question

Get help from others!

© 2024 All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP