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(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^2theta}{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$$$$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")
plt.show()
```

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:

Correct answer by Maxim Umansky on December 9, 2020

Get help from others!

Recent Questions

- How can I transform graph image into a tikzpicture LaTeX code?
- How Do I Get The Ifruit App Off Of Gta 5 / Grand Theft Auto 5
- Iv’e designed a space elevator using a series of lasers. do you know anybody i could submit the designs too that could manufacture the concept and put it to use
- Need help finding a book. Female OP protagonist, magic
- Why is the WWF pending games (“Your turn”) area replaced w/ a column of “Bonus & Reward”gift boxes?

Recent Answers

- Joshua Engel on Why fry rice before boiling?
- haakon.io on Why fry rice before boiling?
- Jon Church on Why fry rice before boiling?
- Peter Machado on Why fry rice before boiling?
- Lex on Does Google Analytics track 404 page responses as valid page views?

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