1

I made a program to calculate the motion of any object (in this case moon) by earth's pull, with zero initial velocity, the moon should just oscillate in a straight line back and forth the earth until it stops at earth exactly right on top of the Earth, so plotting time versus distance I should get a graph similar to a damping signal, instead, I am getting an infinitely decrementing plot.

I have tried :-

  1. giving different initial speeds to the moon, but got the same result, it seems like the way I am using odeint to solve the differential equation is wrong? Not sure. Very new to coding.

  2. Assuming 1000 seconds in not enough for this to happen, so I increased the time to 1e+5,1e+10,1e+20, it seems like odeint couldn't handle that because it gave a different solution every time I run the program for the exact same parameters, received the follow warning:

ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information. warnings.warn(warning_msg, ODEintWarning)

Is there some other function I should use to solve this differential equation?

  1. Reduced the masses to 10,20 and G and r to 10,10, and received the same warning as in case (2)

Any feed back helps

from scipy.integrate import odeint
import matplotlib.pyplot as plt

G=6.67408e-11 #N-m2/kg2 #N-m2/kg2
m1 = 5.972e+24 # kg , mass of earth
m2 = 7.348e+22  # kg , mass of moon


def dvdt(S, t):
    #  v = dr./dt, so a = dv/dt
    r,v = S
    return [v,
            -G*m1 / r**2]

# initial values
r10 = 0 # position of earth in meters
r20 = 4e+8 # position of moon from earth in meters
v10 = 0 # velocity of earth m/s
v20 = 0 # velocity of moon relative to earth m/s
S0 = [r20, v20]
t = np.linspace(0,1000,100)

# solving the differential eqn
acc = odeint(dvdt,S0, t)
r,v = acc.T

# plotting
plt.plot(t,r)
plt.xlabel('Time')
plt.ylabel('Distance between earth and moon')
plt.show()
Mohd. Farhan Hassan
  • 27
  • 2
  • Why do you think a solution would exist past the singularity? This singularity should be reached in finite time, you should get some kind of error message related to the very large acceleration and the correspondingly reducing time step. – Lutz Lehmann Jul 03 '21 at 17:48
  • More of an issue of physics than programming... Does the moon get shot to oblivion when it approaches the earth :D – Futurologist Jul 03 '21 at 19:50
  • @Futurologist haha, not at all, that's the point, that's how I know the plot is wrong. I wrote how it should behave in my question :) – Mohd. Farhan Hassan Jul 04 '21 at 06:03
  • @LutzLehmann, can you explain your comment? It's a simple problem where one body is pulling another one. – Mohd. Farhan Hassan Jul 04 '21 at 06:05
  • You are simulating the free fall straight to the center of gravity of a point mass, with no other forces. The center is reached after a finite time. While the velocity has an upper bound, the acceleration goes to infinity, as do all higher derivatives that are relevant for the step-size control. Try some kind of mollification, the actual gravity force at the center of Earth is zero. – Lutz Lehmann Jul 04 '21 at 06:17
  • @LutzLehmann I tried what you ask, here's what I got ( along with that same warning) image:- https://imgur.com/TdcxuRw – Mohd. Farhan Hassan Jul 04 '21 at 06:47
  • Usually the solver stops at the point of the warning, that is, at radius zero. The remaining data is garbage from a non-zeroed allocation of the result array. // Btw., the force should always point to the center, here that means it should have the opposite sign of the position coordinate. – Lutz Lehmann Jul 04 '21 at 07:27

1 Answers1

0

The scenario assumes that the two bodies are "shifted out of phase", so that they behave like dark matter to each other, no electro-weak or strong nuclear forces.

The effective gravity below the radius of Earth is determined by the mass inside the current radius, the influences of the outer shell add to zero.

The force vector always points to the center, for the one-dimensional motion the sign of the force has always to be opposite to the sign of the coordinate.

In total thus

R1 = 6.7e+6 # m

acc = -sign(r) * G*(m1*min(1,abs(r)/R1)**3) / abs(r)**2
    = -G*m1 * r/max(R1,abs(r))**3

If this is implemented, one gets the expected oscillating plot

plot of fall through the Earth

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
  • @ LutzLehmann Thanks a lot for this, but can you figure out a way using odeint only, the project I am working on requires solving those a lot, so I first took the two body problem as an easy start to slowly progressing into what I really want to do, https://www.youtube.com/watch?v=otRtUiCcCh4&t=706s&ab_channel=Mr.PSolver this is my aim, but I didn't understand much of his math part. – Mohd. Farhan Hassan Jul 04 '21 at 08:16
  • I do not understand your question, I used your code and only changed the acceleration formula as indicated (and switched the grid on). Look up the tag [tag:orbital-mechanics], there are several questions with code for 2-, 3- and N-body systems. – Lutz Lehmann Jul 04 '21 at 08:44