Analytical solution for freefall....

[Will Flannery]

The freefall wiki entry wiki Freefall has an analytic solution for freefall distance in a gravitational field, but ... it doesn't seem to work ... at least i can't get it to work ... here is my MATLAB program to test it ...
clear
G=6.7e-11; % gravitational constant m^3/(kg*s^2)
mEarth = 5.9742e24; % mass of Earth in kg
rEarth = 6.378e6; % radius of Earth in m

% drop a 1kg ball from 100m above the surface of the earth
mu = G*(mEarth + 1);
y0 = rEarth + 100;

t = 1;
x = (3/2*(pi/2-t*sqrt(2*mu/y0^3)))*2/3;
y1 = y0*(x - x^2/5 - 3*x^3/175 - 23*x^4/7875 - 1894*x^5/3931875 - 3293*x^6/21896875 - 2418092*x^7/62077640625);

y0-y1

ans =

9.5639e+04

So, the apple fell 95000 m in the 1st second ... (or I made a mistake)
I tried to check the references, the first has a lot of formulas but not the one above, and the second is behind paywall (The Physics Teacher) and the abstract does not look promising.

Is there an analytic formula anywhere ?

 

[Haborix]

Your definition for x has a typo in the exponent at the end. As an aside, I recall a lot of these series having terrible convergence properties. You might try doing the calculation truncating the inverted series at different powers to see if the answer varies wildly.

 

[Will Flannery]

Thanks, corrected ...
x = (3/2*(pi/2-t*sqrt(2*mu/y0^3)))^2/3;

However, now the apple is falling upwards ...

y0-y1 = ans =

-2.2703e+04

The last term in the series ... 2418092*x^7/62077640625
ans =

0.0028

but ... *y0 = 1.8177e+04

So, it's a long way from converging ... However, as t increases, x decreases and convergence improves, so when t = 500 the last term (including multiplication by y0) is 0.1935 ! So, is it converging to the true value ?

 

[Ibix]

Try ^(2/3). a^2/3 is interpreted as $a^2/3$, but you want $a^{2/3}$.

 

[gmax137]

there are a couple of old threads here where an analytical solution for position vs time is derived. It is messier than one might imagine.

 

[Will Flannery]

Try ^(2/3). a^2/3 is interpreted as $a^2/3$, but you want $a^{2/3}$.

By Jove ! ... and I've only been Matlabbing for 100 years !

G=6.7e-11; % gravitational constant m^3/(kg*s^2)
mEarth = 5.9742e24; % mass of Earth in kg
rEarth = 6.378e6; % radius of Earth in m

mu = G*(mEarth + 1);
y0 = rEarth;

t = 1;
x = (3/2*(pi/2-t*sqrt(2*mu/y0^3)))^(2/3);
y1 = y0*(x - x^2/5 - 3*x^3/175 - 23*x^4/7875 - 1894*x^5/3931875 ...
- 3293*x^6/21896875 - 2418092*x^7/62077640625);

y0-y1 = -2.9134e+04
y0*2418092*x^7/62077640625 = 1.3488e+04

so it's still not converging on the high end ...
for t = 500
y0 -y1 = 1.3217e+06
y0*2418092*x^7/62077640625 = 296.6925

So, it's better but still not so good.
Conclusion - there are not enough terms in this thing to tell if it works or not (or there is yet another bug ... )

 

[vanhees71]

Just use LaTeX, and one can read you math!

https://www.physicsforums.com/help/latexhelp/

 

[Will Flannery]

Just use LaTeX, and one can read you math!

https://www.physicsforums.com/help/latexhelp/

That's not math, that's MATLAB.

In any case, now the question is, why does every book ignore the problem of an analytical solution for the 1-D gravity problem (as noted in the ref. in the wiki article, the guy looked at 100 physics books) and what is the solution?

From the reference in the wiki article, the approach to a solution is thru Kepler's eq. Note that when eccentricity e approaches 1 in Kepler's eq, the orbit does not approach a parabola, Kepler's original eq. only works for ellipses, and as e -> 1 the orbit gets flatter and flatter and when e = 1 you get a straight line !, I just plotted one, So, apparently we can use either a numeric or analytic solution to Kepler's eq, with e = 1, to plot a 1-d trajectory. And so the questions for the day is ... how?

So, if we drop a ball from a height of 100m, a = 100+rEarth and we can plot the trajectory ...
(1.5 hours later) ... that seems to work ...
clear
G=6.7e-11; % gravitational constant m^3/(kg*s^2)
mEarth = 5.9742e24; % mass of Earth in kg
rEarth = 6.378e6; % radius of Earth in m

a = 100+rEarth;
T = sqrt(a^3 /(G*mEarth/(4*pi^2))); % from Kepler's 3rd Law
e = 1;

N=100;
for i = 1:4
E = 2*pi*(i-1)/N + pi; % Eccentric anomaly, start at apogee
M = E - e*sin(E); % Mean anomaly
t(i) = M*T/(2*pi)-T/2; % M(t) = 2*pi*t/T
x(i) = a*cos(E);
end

plot(t,-x-rEarth) % invert graph

 

[vanhees71]

Why are you writing when you don't want to be read? SCNR.

 

[Will Flannery]

Why are you writing when you don't want to be read? SCNR.

? That doesn't make sense, I want it to be read... even debugged if necessary ...

Generally, things become much clearer when you actually have to program them and check if they work ... For example, it took me over an hour, and many corrections, to get that program to work.

 

[Ibix]

If it's code, put it in code tags.

X=1 % It's even got a syntax highlighter

 

[gmax137]

If it's code, put it in code tags.

SCNR

 

[jasonRF]

EDIT: after thinking for a few more minutes, I think my concern is bogus. I think my post isn't worth reading, but don't like to delete posts since it feels dishonest...

I know this is a very old thread, but I just stumbled across it. If you still care about this it seems to me that you have a precision problem. Matlab uses double precision by default, so mEarth+1=mEarth. The smallest mass you could use and not get this problem will be the order of 1e9 kg. Now the value of $\mu$ doesn't change the result, but you may have similar problems with other portions of this calculation as well. Have you looked at the orders of magnitude of each of the terms in your sum?

I have certainly run into cases like this in my work on numerous occasions, where a pretty analytical formula is completely useless for practical calculations. Sometimes using appropriate numerical techniques is much more practical.


jason