Determining effect of small variable force on planetary perihelion precession Is there an analytical technique for determining the effect of a small variable transverse acceleration upon the rate of aspides precession (strictly not a precession but rotation of the line of aspides) of a planet orbitting around the Sun in a 2D plane according to Newtonian gravity law?  

I have modelled such effects in a reiterative computer model and would like to verify those measurements.  

The transverse acceleration formula is 


  $$At = (K/c^2)*Vr*Vt * Ar.$$


Where:-


  c is speed of light, 
  
  K is a constant of magnitude between 0 and +/-3, such that $K/(c^2) &lt;&lt; 1$. 
  
  Ar is the acceleration of the planet towards the Sun due to Newtonian gravitational influence of the Sun, ($Ar = GM/r^2$).
  
  Vr is radial component of planet velocity relative to the Sun (+ = motion away from the Sun)
  
  Vt is transverse  component of planet velocity relative to the Sun (+ = direction of planet forward motion along its orbital path). Vectorially Vt = V - Vr where V is the total instantaneous velocity vector of the planet relative to the Sun. 
  
  Assume planet mass is small relative to the Sun
  
  No other bodies are in the system
  
  All motions and accelerations are confined to the two-dimensional plane of the orbit.


UPDATE

The reason why this is interesting to me is that a value of K = +3 in my computer model produces anomalous (Non-Newtonian) periapse rotation rate values very close within about 1% of those predicted by General Relativity and within a few percent of those observed by astronomers (Le Verrier, updated by Newcomb).

Formula (Einstein, 1915) for GR-derived periapse rotation (radians per orbit)
from http://en.wikipedia.org/wiki/Apsidal_precession
$$
\boldsymbol{\omega}=24.\pi^3.a^2.T^{-2}.c^{-2}.(1-e^2)^{-1}
$$

UPDATE 4

I have accepted Walter's answer.  Not only did he answer the original question (Is there a technique...?) but also his analysis produces a formula which not only confirms the computer-simulated effects of the transverse acceleration formula (for K = 3) but which also (unexpectedly to me) is essentially equivalent to the Einstein 1915 formula.

from Walter's Summary (in Walter's answer below):-

: (from first order peturbation analysis) semi-major axis and eccentricity are unchanged, but the direction of periapse rotates in the plane of the orbit at rate $$ \omega=\Omega \frac{v_c^2}{c^2} \frac{K}{1-e^2}, $$ where $\Omega$ is the orbital frequency and $v_c=\Omega a$ with $a$ the semi-major axis. Note that (for $K=3$) this agrees with the general relativity (GR) precession rate at order $v_c^2/c^2$ (given by Einstein 1915).
You may want to use perturbation theory. This only gives you an approximate answer, but allows for analytic treatment. Your force is considered a small perturbation to the Keplerian elliptic orbit and the resulting equations of motion are expanded in powers of $K$. For linear perturbation theory, only terms linear in $K$ are retained. This simply leads to integrating the perturbation along the unperturbed original orbit. Writing your force as a vector, the perturbing acceleration is
$$
\boldsymbol{a} = K \frac{GM}{r^2c^2}v_r\boldsymbol{v}_t
$$
with $v_r=\boldsymbol{v}{\cdot}\hat{\boldsymbol{r}}$ the radial velocity 
($\boldsymbol{v}\equiv\dot{\boldsymbol{r}}$) and
$\boldsymbol{v}_t=(\boldsymbol{v}-\hat{\boldsymbol{r}}(\boldsymbol{v}{\cdot}\hat{\boldsymbol{r}}))$ the rotational component of velocity (the full velocity minus the radial velocity). Here, the dot above denotes a time derivative and a hat the unit vector.

Now, it depends what you mean with 'effect'. Let's work out the changes of the orbital semimajor axis $a$, eccentricity $e$, and direction of periapse. 



To summarise the results below: semi-major axis and eccentricity are unchanged, but the direction of periapse rotates in the plane of the orbit at rate
$$
\omega=\Omega \frac{v_c^2}{c^2}
\frac{K}{1-e^2},
$$
where $\Omega$ is the orbital frequency and $v_c=\Omega a$ with $a$ the semi-major axis. Note that (for $K=3$) this agrees with the general relativity (GR) precession rate at order $v_c^2/c^2$ (given by Einstein 1915 but not mentioned in the original question).



change of semimajor axis

From the relation $a=-GM/2E$ (with $E=\frac{1}{2}\boldsymbol{v}^2-GMr^{-1}$ the orbital energy) we have for the change of $a$ due to an external (non-Keplerian) acceleration
$$
\dot{a}=\frac{2a^2}{GM}\boldsymbol{v}{\cdot}\boldsymbol{a}.
$$
Inserting $\boldsymbol{a}$ (note that $\boldsymbol{v}{\cdot}\boldsymbol{v}_t=h^2/r^2$ with angular momentum vector $\boldsymbol{h}\equiv\boldsymbol{r}\wedge\boldsymbol{v}$), we get
$$
\dot{a}=\frac{2a^2Kh^2}{c^2}\frac{v_r}{r^4}.
$$
Since the orbit average $\langle v_r f(r)\rangle=0$ for any function $f$ (see below), $\langle\dot{a}\rangle=0$. 

change of eccentricity

From $\boldsymbol{h}^2=(1-e^2)GMa$, we find
$$
e\dot{e}=-\frac{\boldsymbol{h}{\cdot}\dot{\boldsymbol{h}}}{GMa}+\frac{h^2\dot{a}}{2GMa^2}.
$$
We already know that $\langle\dot{a}\rangle=0$, so only need to consider the first term. Thus,
$$
e\dot{e}=-\frac{(\boldsymbol{r}\wedge\boldsymbol{v}){\cdot}(\boldsymbol{r}\wedge\boldsymbol{a})}{GMa}
         =-\frac{r^2\;\boldsymbol{v}{\cdot}\boldsymbol{a}}{GMa}
         =-\frac{Kh^2}{ac^2}\frac{v_r}{r^2},
$$
where I have used the identity
$(\boldsymbol{a}\wedge\boldsymbol{b}){\cdot}(\boldsymbol{c}\wedge\boldsymbol{d})
=\boldsymbol{a}{\cdot}\boldsymbol{c}\;\boldsymbol{b}{\cdot}\boldsymbol{d}-
 \boldsymbol{a}{\cdot}\boldsymbol{d}\;\boldsymbol{b}{\cdot}\boldsymbol{c}$ and the fact $\boldsymbol{r}{\cdot}\boldsymbol{a}_p=0$.
Again $\langle v_r/r^2\rangle=0$ and hence $\langle\dot{e}\rangle=0$.

change of the direction of periapse

The eccentricity vector
$
\boldsymbol{e}\equiv\boldsymbol{v}\wedge\boldsymbol{h}/GM - \hat{\boldsymbol{r}}
$
points (from the centre of gravity) in the direction of periapse, has magnitude $e$, and is conserved under the Keplerian motion (validate all that as an exercise!). From this definition we find its instantaneous change due to external acceleration
$$
\dot{\boldsymbol{e}}=
\frac{\boldsymbol{a}\wedge(\boldsymbol{r}\wedge\boldsymbol{v})
     +\boldsymbol{v}\wedge(\boldsymbol{r}\wedge\boldsymbol{a})}{GM}
     =\frac{2(\boldsymbol{v}{\cdot}\boldsymbol{a})\boldsymbol{r}
       -(\boldsymbol{r}{\cdot}\boldsymbol{v})\boldsymbol{a}}{GM}
     =\frac{2K}{c^2}\frac{h^2v_r\boldsymbol{r}}{r^4}
     -\frac{K}{c^2}\frac{v_r^2\boldsymbol{v}_t}{r}
$$
where I have used the identity 
$\boldsymbol{a}\wedge(\boldsymbol{b}\wedge\boldsymbol{c})=(\boldsymbol{a}{\cdot}\boldsymbol{c})\boldsymbol{b}-(\boldsymbol{a}{\cdot}\boldsymbol{b})\boldsymbol{c}$
and the fact $\boldsymbol{r}{\cdot}\boldsymbol{a}=0$. The orbit averages of these expression are considered in the appendix below. If we finally put everything together, we get
$
\dot{\boldsymbol{e}}=\boldsymbol{\omega}\wedge\boldsymbol{e}
$
with [corrected again]
$$
\boldsymbol{\omega}=\Omega K \frac{v_c^2}{c^2}
(1-e^2)^{-1}\, \hat{\boldsymbol{h}}.
$$
This is a rotation of periapse in the plane of the orbit with angular frequency $\omega=|\boldsymbol{\omega}|$. In particular $\langle e\dot{e}\rangle=\langle\boldsymbol{e}{\cdot}\dot{\boldsymbol{e}}\rangle=0$ in agreement with our previous finding.

Don't forget that due to our usage of first-order perturbation theory these results are only strictly true in the limit $K(v_c/c)^2\to0$. At second-order perturbation theory, however, both $a$ and/or $e$ may change. In your numerical experiments, you should find that the orbit-averaged changes of $a$ and $e$ are either zero or scale stronger than linear with perturbation amplitude $K$.

disclaimer No guarantee that the algebra is correct. Check it!



Appendix: orbit averages

Orbit averages of $v_rf(r)$ with an abitrary (but integrable) function $f(r)$ can be
directly calculated for any type of periodic orbit. Let $F(r)$ be the antiderivative of $f(r)$, i.e. $F'\!=f$, then the orbit average is:
$$
\langle v_r f(r)\rangle = \frac{1}{T}\int_0^T v_r(t)\,f\!\left(r(t)\right) \mathrm{d}t
= \frac{1}{T} \left[F\left(r(t)\right)\right]_0^T = 0
$$
with $T$ the orbital period.

For the orbit averages required in $\langle\dot{\boldsymbol{e}}\rangle$, we must dig a bit deeper. For a Keplerian elliptic orbit
$$
\boldsymbol{r}=a\left((\cos\eta-e)\hat{\boldsymbol{e}}+\sqrt{1-e^2}\sin\eta\,\hat{\boldsymbol{k}}\right)\qquad\text{and}\qquad
r=a(1-e\cos\eta)
$$
with eccentricity vector $\boldsymbol{e}$ and $\hat{\boldsymbol{k}}\equiv\hat{\boldsymbol{h}}\wedge\hat{\boldsymbol{e}}$ a vector perpendicular to $\boldsymbol{e}$ and $\boldsymbol{h}$. Here, $\eta$ is the eccentric anomaly, which is related to the mean anomaly $\ell$ via
$
\ell=\eta-e\sin\eta,
$
such that $\mathrm{d}\ell=(1-e\cos\eta)\mathrm{d}\eta$ and an orbit average becomes
$$
\langle\cdot\rangle = (2\pi)^{-1}\int_0^{2\pi}\cdot\;\mathrm{d}\ell = (2\pi)^{-1}\int_0^{2\pi}\cdot\;(1-e\cos\eta)\mathrm{d}\eta.
$$
Taking the time derivative (note that $\dot{\ell}=\Omega=\sqrt{GM/a^3}$ the orbital frequency) of $\boldsymbol{r}$, we find for the instantaneous (unperturbed) orbital velocity
$$
\boldsymbol{v}=v_c\frac{\sqrt{1-e^2}\cos\eta\,\hat{\boldsymbol{k}}-\sin\eta\,\hat{\boldsymbol{e}}}{1-e\cos\eta}
$$
where I have introduced $v_c\equiv\Omega a=\sqrt{GM/a}$, the speed of the circular orbit with semimajor axis $a$. From this, we find the radial velocity $v_r=\hat{\boldsymbol{r}}{\cdot}\boldsymbol{v}=v_c e\sin\eta(1-e\cos\eta)^{-1}$
and the rotational velocity
$$
\boldsymbol{v}_t = v_c\frac{\sqrt{1-e^2}(\cos\eta-e)\,\hat{\boldsymbol{k}}-(1-e^2)\sin\eta\,\hat{\boldsymbol{e}}}{(1-e\cos\eta)^2}.
$$

With these, we have [corrected again]
$$
\left\langle \frac{h^2v_r\boldsymbol{r}}{r^4}\right\rangle = 
\Omega v_c^2\,\hat{\boldsymbol{k}}\,
\frac{e(1-e^2)^{3/2}}{2\pi}\int_0^{2\pi}\frac{\sin^2\!\eta}{(1-e\cos\eta)^4}\mathrm{d}\eta
=\frac{\Omega v_c^2e}{2(1-e^2)}\hat{\boldsymbol{k}}
\\
\left\langle \frac{v_r^2\boldsymbol{v}_t}{r}\right\rangle = \Omega v_c^2\,
\hat{\boldsymbol{k}}\,
\frac{e^2(1-e^2)^{1/2}}{2\pi}\int_0^{2\pi}\frac{\sin^2\!\eta(\cos\eta-e)}{(1-e\cos\eta)^4}\mathrm{d}\eta=0,
$$
in particular, the components in direction $\hat{\boldsymbol{e}}$ average to zero. Thus [corrected again]
$$\left\langle 2\frac{h^2v_r\boldsymbol{r}}{r^4}-\frac{v_r^2\boldsymbol{v}_t}{r}\right\rangle
=\frac{\Omega v_c^2e\,\hat{\boldsymbol{k}}}{(1-e^2)}
$$