Solving Kepler's 1st law as a function of time

[jebez]

TL;DR Summary

p , ε , c constants :

θ'(t)(p/(1+ε cos(θ(t))))^2=c

θ(t)=?

Hi

I posted this differential equation to WolframAlpha https://www.wolframalpha.com/input?i2d=true&i=Power[\(40)Divide[a,1+b*cos\(40)y\(40)x\(41)\(41)]\(41),2]*y'\(40)x\(41)=c but no solution , " Standard computation time exceeded... Try again with Pro computation time "
Should I ( buy and ) post to Wolfram|Alpha Pro ? I don't want to buy Pro if it doesn't solve it ... Maybe someone who has Pro can post it ? Lend his/her account ?

In fact we have r(θ) = p/(1+ε cos(θ)) ( 1st Kepler's law ) , r²θ'(t)=c ( 2nd Kepler's law ) and I want r(t) so need to solve the differential equation .
https://en.wikipedia.org/wiki/Kepler's_laws_of_planetary_motion

Thanks .

 

[PeroK]

In fact we have r(θ) = p/{1 + ε cos(θ)) ( 1st Kepler law ) , r² θ'(t) = c ( 2nd Kepler law ) and I want r(t) so need to solve the differential equation .
https://en.wikipedia.org/wiki/Kepler's_laws_of_planetary_motion

Thanks .

Surprisingly, perhaps, there is no closed-form solution other than for circular orbits.

 

[fresh_42]

TL;DR Summary: p , ε , c constants :

θ'(t) (p/(1 + ε cos(θ(t))))^2) = c

θ(t) = ?

Hi

I posted this differential equation to WolframAlpha https://www.wolframalpha.com/input?i2d=true&i=Power[\(40)Divide[a,1+b*cos\(40)y\(40)x\(41)\(41)]\(41),2]*y'\(40)x\(41)=c but no solution , " Standard computation time exceeded... Try again with Pro computation time "
Should I ( buy and ) post to Wolfram|Alpha Pro ? I don't want to buy Pro if it doesn't solve it ... Maybe someone who has Pro can post it ? Lend his/her account ?

In fact we have r(θ) = p/(1 + ε cos(θ)) ( 1st Kepler law ) , r² θ'(t) = c ( 2nd Kepler law ) and I want r(t) so need to solve the differential equation .
https://en.wikipedia.org/wiki/Kepler's_laws_of_planetary_motion

Thanks .

Have a look at
https://www.wolframalpha.com/input?i=arccos(y(t))'=(a+b*+y(t))^2

 

[pasmith]

The equation is separable in $\theta$, so $$\int_{\theta(0)}^{\theta(t)} \frac{1}{(1 + \epsilon\cos\phi)^2}\,d\phi = \frac{ct}{p^2}.$$ But that's as far as you can go: there is no known antiderivative for the integrand when $\epsilon \neq 0$.

 

[fresh_42]

The equation is separable in $\theta$, so $$\int_{\theta(0)}^{\theta(t)} \frac{1}{(1 + \epsilon\cos\phi)^2}\,d\phi = \frac{ct}{p^2}.$$ But that's as far as you can go: there is no known antiderivative for the integrand when $\epsilon \neq 0$.

What's wrong with my substitution $y(t)=\cos\theta(t)$? Nasty, but it looks like it has an implicit solution.

 

[jebez]

Well sorry I thought equations will be rendered automatically like WolframAlpha , it seems to need LaTeX but I never used this and it seems tiresome ...

In fact I want to verify if r(t) is solution to this differential equation to prove the Kepler's 1st law :

c , G , M constants :

c²/r(t)^3-G M/r(t)²-r''(t)=0

r(t)=?
r(θ)=?

Having again r²θ'(t)=c ( 2nd Kepler's law ) , centrifugal force=m r θ'(t)² and gravity force=G M m/r² so
force=m r''(t)=centrifugal force - gravity force
m r''(t)=m c²/r(t)^3-G M m/r²
then above .

And same problem with WolframAlpha .

 

[pasmith]

The usual approach is to set $u = 1/r$ and $$\begin{split} \frac{dr}{dt} &= \frac{dr}{d\theta}\frac{d\theta}{dt} = \frac{L}{r^2}\frac{dr}{d\theta} = -L\frac{du}{d\theta} \\ \frac{d^2r}{dt^2} &= \frac{L}{r^2} \frac{d}{d\theta}\frac{dr}{dt} = -L^2u^2\frac{d^2 u}{d\theta^2} \end{split}$$ whence $$\frac{d^2u}{d\theta^2} + u = \frac{GM}{L^2}$$ which we can solve for $u(\theta)$.

 

[pasmith]

What's wrong with my substitution $y(t)=\cos\theta(t)$? Nasty, but it looks like it has an implicit solution.

In fact WolframAlpha does give an antiderivative $$\int \frac{1}{(1 + \epsilon\cos\theta)^2}\,d\theta = \frac{\epsilon \sin \theta}{(\epsilon^2 - 1)(\epsilon \cos \theta + 1)} - \frac{2\operatorname{artanh}\left( \frac{(\epsilon-1)\tan(\frac12\theta)}{\sqrt{\epsilon^2-1}} \right)}{(\epsilon^2 - 1)^{3/2}} + C$$ but we then have the problem of solving that for $\theta$.

 

[jebez]

I found r=p/(1+ε cos(θ(t)) is solution to c²/r(t)^3-G M/r(t)²-r''(t)=0 :

r²θ'(t)=c
θ'(t)=c/r²
(dr/dθ)(dθ/dt)=dr/dt=v=r'(θ)c/r²=c ε sin(θ)/p

v'(θ)=c ε cos(θ)/p
(dv/dθ)(dθ/dt)=dv/dt=r''(t)=ε cos(θ)(c(/1+ε cos(θ)))²/p^3

c²/r(t)^3-G M/r(t)²-r''(t)=0
c²(1+ε cos(θ)/p)^3-G M(1+ε cos(θ)/p)²-ε cos(θ)(c(/1+ε cos(θ)))²/p^3=0

c²-p G M=0 as bonus .

Feel free to edit my posts to render equations , yes I can do it myself but it would be nice automatically ...

 

[berkeman]

Feel free to edit my posts to render equations , yes I can do it myself but it would be nice automatically ...

Seriously?

Which word in the "LaTeX Guide" link below the Edit window did you have trouble understanding?

 

[jebez]

I found t(θ) :
1/θ'(t)=dt/dθ=t'(θ)=(p/(1+ε cos(θ)))²/c
t(θ)=see https://www.wolframalpha.com/input?i2d=true&i=Divide[a²*antiderivative+\(40)Divide[1,1+b*cos\(40)x\(41)]\(41)²,c] .

Now we need to from t(θ) to θ(t) , how to do this ? With WolframAlpha ?

I can't edit my previous messages to render equations , I'm just lazy to do that sorry and it isn't copyable ...

 

[berkeman]

I can't edit my previous messages to render equations , I'm just lazy to do that sorry and it isn't copyable ...

Please be sure to use LaTeX from now on here. Thank you.

 

[pasmith]

- Copy your

Well sorry I thought equations will be rendered automatically like WolframAlpha , it seems to need LaTeX but I never used this and it seems tiresome ...

In fact I want to verify if r(t) is solution to this differential equation to prove the Kepler's 1st law :

c , G , M constants :

c²/r(t)^3-G M/r(t)²-r''(t)=0

As I suggested above, you should substitute $u = 1/r$ and eliminate $t$ as the independent variable in favour of $\theta$. This gives you $$u'' + u = \frac{GM}{L^2}$$ where $L = r^2 \dot \theta$ is constant. This is a linear ODE with constant coefficients which can be easily solved: $$r(\theta) =\frac{1}{u(\theta)} = \frac{L^2/GM}{1 + \epsilon\cos(\theta - \theta_0)}.$$

 

[jebez]

I can't edit my messages #6 and #9 but r=c²/(G M(1+ε cos(θ))) , c²/r(t)³-G M/r(t)²-r''(t)=0 can be simplified :
c²/r(t)-r(t)²r''(t)=G M
c²G M(1+ε cos(θ))/c²-c⁴/(G M(1+ε cos(θ)))²ε cos(θ)(G M)³(1+ε cos(θ))²/c⁴=G M .

Sorry pasmith but I don't understand from -L du/dθ in #7 ...

Well I found how to inverse t(θ) to θ(t) with WolframAlpha and Mathematica :

InverseFunction[(a^2 Integrate[(1/(1 + b Cos[x]))^2, x])/c]

but both didn't solve it , see https://community.wolfram.com/groups/-/m/t/2697668?p_p_auth=2VkMYTxA .

And ( maybe ) simplification :

t(θ)=(a^2 ((b sin(x))/((b^2-1) (b cos(x)+1))-(2 (b-1) tan(x/2) coth(x))/(b^2-1)^2+d))/c .

 

[jebez]

t(θ)=$$\frac{a^2 \left(\frac{b (x \sin )}{\left(b^2-1\right) (b (x \cos )+1)}-\frac{2 (b-1) (x \tan ) (x \coth )}{2 \left(b^2-1\right)^2}+d\right)}{c}$$
( a=p , b=ε , c=c , d=t₀ )
given by Mathematica ( I've 15 days trial ) , not WolframAlpha .

There's also the differential equation θ'(t)=c((1+ε cos(θ(t)))/p)² aka y'(x)=c((1+b cos(y(x)))/a)² in Mathematica gives :

{{y(x)->InverseFunction[1/2 ((b sin(#1))/((b^2-1) (b cos(#1)+1))-(2 tanh^-1(((b-1) tan(#1/2))/Sqrt[b^2-1]))/(b^2-1)^(3/2))&][(c x)/(2 a^2)+Subscript[c, 1]]}}

$$\left\{\left\{y(x)\to \text{InverseFunction}\left[\frac{1}{2} \left(\frac{b (\text{$\#$1} \sin )}{\left(b^2-1\right) (b (\text{$\#$1} \cos )+1)}-\frac{2 ((b-1) (\text{$\#$1} \tan ))}{\left(b^2-1\right)^{3/2} \left(\left(2 \sqrt{b^2-1}\right) \tanh \right)}\right)\&\right]\left[\frac{c x}{2 a^2}+c_1\right]\right\}\right\}$$
What is # and & here ?

With Mathematica I can copy as LaTeX .

 

[fresh_42]

With Mathematica I can copy as LaTeX .

https://www.physicsforums.com/help/latexhelp/
http://detexify.kirelabs.org/symbols.html

are basically all you need.

 

[jebez]

I'm back , sorry for long time .

In short :

$L=r^2ω=constant$ ( Kepler's 2nd law )

$r=\frac{L^2}{GM(1+εcos(θ))}$ ( Kepler's 1st law )

⇒$\frac{1}{ω}=\frac{dt}{dθ}=\frac{L^3}{(GM(1+εcos(θ)))^2}$

⇔$t(θ)=\frac{L^3}{(GM)^2}∫\frac{1}{(1+εcos(θ))^2}dθ$


$∫\frac{1}{(1+εcos(θ))^2}dθ$ :

SageMath 10.3 via CoCalc :

In :
var('θ,ε')
assume(0<ε<1)
from sage.symbolic.integration.integral import indefinite_integral
latex(indefinite_integral(1/(1+ε*cos(θ))^2,θ))

Out :
$\frac{2 \, ε \sin\left(θ\right)}{{\left(ε^{3} + ε^{2} - ε - \frac{{\left(ε^{3} - ε^{2} - ε + 1\right)} \sin\left(θ\right)^{2}}{{\left(\cos\left(θ\right) + 1\right)}^{2}} - 1\right)} {\left(\cos\left(θ\right) + 1\right)}} + \frac{2 \, \arctan\left(\frac{{\left(ε - 1\right)} \sin\left(θ\right)}{\sqrt{-ε^{2} + 1} {\left(\cos\left(θ\right) + 1\right)}}\right)}{{\left(ε^{2} - 1\right)} \sqrt{-ε^{2} + 1}}$


$t(θ)=\frac{2L^3}{(GM)^2}\left(\frac{ε \sin\left(θ\right)}{{\left(ε^{3} + ε^{2} - ε - \frac{{\left(ε^{3} - ε^{2} - ε + 1\right)} \sin\left(θ\right)^{2}}{{\left(\cos\left(θ\right) + 1\right)}^{2}} - 1\right)} {\left(\cos\left(θ\right) + 1\right)}} + \frac{\arctan\left(\frac{{\left(ε - 1\right)} \sin\left(θ\right)}{\sqrt{-ε^{2} + 1} {\left(\cos\left(θ\right) + 1\right)}}\right)}{{\left(ε^{2} - 1\right)} \sqrt{-ε^{2} + 1}}\right)+t_0$


For θ(t) , inversing a function in SageMath isn't yet implemented ...
https://ask.sagemath.org/question/8071/can-sage-compute-the-inverse-of-a-function/
https://github.com/sagemath/sage/issues/11202

 

[jebez]

I can't edit my previous post ...

There're (cos(θ)+1) as denominators , so we must precise :

-π<θ<π

 

[jebez]

In :
var('θ,ε,t')
assume(0<ε<1)
assume(-π<θ<π)
from sage.symbolic.integration.integral import indefinite_integral
latex(solve(t==indefinite_integral(1/(1+ε*cos(θ))^2,θ),θ))

Out :
$\left[\sin\left(θ\right) = \frac{\sqrt{-ε^{2} + 1} {\left(ε \cos\left(θ\right) + ε\right)} - \sqrt{-t^{2} ε^{8} + 4 \, t^{2} ε^{6} - {\left(6 \, t^{2} + 1\right)} ε^{4} + {\left(4 \, t^{2} + 1\right)} ε^{2} + 4 \, {\left(ε^{2} - 1\right)} \arctan\left(\frac{{\left(ε - 1\right)} \sin\left(θ\right)}{\sqrt{-ε^{2} + 1} {\left(\cos\left(θ\right) + 1\right)}}\right)^{2} - 4 \, {\left(t ε^{4} - 2 \, t ε^{2} + t\right)} \sqrt{-ε^{2} + 1} \arctan\left(\frac{{\left(ε - 1\right)} \sin\left(θ\right)}{\sqrt{-ε^{2} + 1} {\left(\cos\left(θ\right) + 1\right)}}\right) - t^{2}} {\left(\cos\left(θ\right) + 1\right)}}{2 \, {\left(ε - 1\right)} \arctan\left(\frac{{\left(ε - 1\right)} \sin\left(θ\right)}{\sqrt{-ε^{2} + 1} {\left(\cos\left(θ\right) + 1\right)}}\right) - {\left(t ε^{3} - t ε^{2} - t ε + t\right)} \sqrt{-ε^{2} + 1}}, \sin\left(θ\right) = \frac{\sqrt{-ε^{2} + 1} {\left(ε \cos\left(θ\right) + ε\right)} + \sqrt{-t^{2} ε^{8} + 4 \, t^{2} ε^{6} - {\left(6 \, t^{2} + 1\right)} ε^{4} + {\left(4 \, t^{2} + 1\right)} ε^{2} + 4 \, {\left(ε^{2} - 1\right)} \arctan\left(\frac{{\left(ε - 1\right)} \sin\left(θ\right)}{\sqrt{-ε^{2} + 1} {\left(\cos\left(θ\right) + 1\right)}}\right)^{2} - 4 \, {\left(t ε^{4} - 2 \, t ε^{2} + t\right)} \sqrt{-ε^{2} + 1} \arctan\left(\frac{{\left(ε - 1\right)} \sin\left(θ\right)}{\sqrt{-ε^{2} + 1} {\left(\cos\left(θ\right) + 1\right)}}\right) - t^{2}} {\left(\cos\left(θ\right) + 1\right)}}{2 \, {\left(ε - 1\right)} \arctan\left(\frac{{\left(ε - 1\right)} \sin\left(θ\right)}{\sqrt{-ε^{2} + 1} {\left(\cos\left(θ\right) + 1\right)}}\right) - {\left(t ε^{3} - t ε^{2} - t ε + t\right)} \sqrt{-ε^{2} + 1}}\right]$

t(θ) seems non-invertible .

 

[Will Flannery]

From the title, I think you want to be able to calculate orbit position as a function of time?

That's known as Kepler's Problem, which is to solve Kepler's Equation, which can be derived from Kepler's 3 laws.

Kepler's equation - https://en.wikipedia.org/wiki/Kepler's_equation

Kepler's problem - https://en.wikipedia.org/wiki/Kepler_problem
The force may be either attractive or repulsive. The problem is to find the position or speed of the two bodies over time given .......

Funny story - You'll notice that the wiki entry for the Kepler problem has a section Solution of the Kepler Problem, but the equation derived is not a solution to the Kepler problem, as it calculates u = 1/radius as a function of an angle.

The reference for Kepler's problem is the book - https://www.amazon.com/Solving-Keplers-Equation-Three-Centuries/dp/0943396409?tag=pfamazon01-20

 

[jebez]

https://en.wikipedia.org/wiki/Kepler's_equation just gives at the end the Kepler's 1st law of my #17 post in general ( $-\frac{k}{m}$instead of $GM$ ) .

The problem is solved , we've r(θ) & t(θ) , the summary with " the conversation ultimately ends without a solution " needs to be edited .

Now why not put it on Wikipedia ? Any volunteer ?

https://www.amazon.com/Solving-Keplers-Equation-Three-Centuries/dp/0943396409?tag=pfamazon01-20 : $149 more than Wikipedia but well I'll buy it , thanks for the link ( irony , sarcasm ) .

 

[Will Flannery]

jebez: (irony, sarcasm) ...... I didn't notice the price! I bought the book several years ago and it wasn't expensive, I just checked and found 1 copy ! for #30.00 https://shopatsky.com/products/solving-keplers-equation

Inverting Kepler's equation (not Kepler's 1st law), i.e. calculating r(t), is known as Kepler's problem. There are two ways to solve it, numerically and analytically.

Both Kepler and Newton solved it numerically, Kepler's method was limited to not-too-eccentric orbits, Newton's method is general. The wiki entry for Kepler's equation contains a section on Newton's method.

Solving it analytically is a whole new ballgame. Lagrange derived an analytic solution in the 1700s, and it's been a hot topic in some circles since that time, i.e. three centuries. The analytic solutions are infinite series and converge for some values and don't converge for others, i.e. they can be problematical. This is the subject of Colwell's book. It was a revelation to me (and I had worked several years analyzing orbits and rocket trajectories in the aerospace industry).

The only place I've seen an analytic solution online is https://en.wikipedia.org/wiki/Free_fall which has a 1-d solution in the section on an 'inverse square-law gravitational field.'.

 

[jebez]

Erratum :
https://en.wikipedia.org/wiki/Kepler_problem instead of https://en.wikipedia.org/wiki/Kepler's_equation in my #21 post .

In Maxima :

assume(ε>0,ε<1);
integrate(1/(1+ε*cos(θ))^2,θ);

assume(ε>0,ε<1,θ>-π,θ<π);
solve(t=integrate(1/(1+ε*cos(θ))^2,θ),θ);

 

[jebez]

Earth's elliptic orbit :

assume(ε>0,ε<1);
plot2d(2456084229048821270250000000/489238659320312440009*subst(ε=0.0167086,integrate(1/(1+ε*cos(θ))^2,θ)),[θ,-π,π]);