<p>You may want to use <strong>perturbation theory</strong>. This only gives you an <em>approximate</em> 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.</p>

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

<hr>

<p>To <strong>summarise the results below</strong>: 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 <strong>agrees with the general relativity (GR) precession rate</strong> at order $v_c^2/c^2$ (given by Einstein 1915 but not mentioned in the original question).</p>

<hr>

<p><strong>change of semimajor axis</strong></p>

<p>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$. </p>

<p><strong>change of eccentricity</strong></p>

<p>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$.</p>

<p><strong>change of the direction of periapse</strong></p>

<p>The <em>eccentricity vector</em>
$
\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 [<em>corrected again</em>]
$$
\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.</p>

<p>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$.</p>

<p><strong>disclaimer</strong> No guarantee that the algebra is correct. Check it!</p>

<hr>

<p><strong>Appendix: orbit averages</strong></p>

<p>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.</p>

<p>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}.
$$</p>

<p>With these, we have [<em>corrected again</em>]
$$
\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 [<em>corrected again</em>]
$$\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)}
$$</p>

CC BY-SA 3.0
2014-09-04T21:16:02.710Z
Comments are not for extended discussion; this conversation has been [moved to chat](http://chat.stackexchange.com/rooms/16932/discussion-on-answer-by-walter-determining-effect-of-small-variable-force-on-pla).
called2voyage
CC BY-SA 3.0
2014-09-01T19:27:42.863Z
2014-09-09T08:35:49.607Z
2014-09-09T08:35:49.607Z
Walter
<p>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?  </p>

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

<p>The transverse acceleration formula is </p>

<blockquote>
  <p>$$At = (K/c^2)*Vr*Vt * Ar.$$</p>
</blockquote>

<p>Where:-</p>

<blockquote>
  <p><em>c</em> is speed of light, </p>
  
  <p>K is a constant of magnitude between 0 and +/-3, such that $K/(c^2) &lt;&lt; 1$. </p>
  
  <p>Ar is the acceleration of the planet towards the Sun due to Newtonian gravitational influence of the Sun, ($Ar = GM/r^2$).</p>
  
  <p>Vr is radial component of planet velocity relative to the Sun (+ = motion away from the Sun)</p>
  
  <p>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. </p>
  
  <p>Assume planet mass is small relative to the Sun</p>
  
  <p>No other bodies are in the system</p>
  
  <p>All motions and accelerations are confined to the two-dimensional plane of the orbit.</p>
</blockquote>

<p>UPDATE</p>

<p>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).</p>

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

<p>UPDATE 4</p>

<p>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.</p>

<p>from Walter's Summary (in Walter's answer below):-</p>

<p>: (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).</p>

CC BY-SA 3.0
2014-08-27T15:52:52.647Z
Are you still seeking an answer?
Walter
CC BY-SA 3.0
2014-09-01T00:10:46.790Z
@Walter.  Yes I am.  I have asked similar question at https://physics.stackexchange.com/questions/123685/effect-of-small-variable-transverse-force-on-satellite-orbital-precession but no solid answer received yet.
steveOw
CC BY-SA 3.0
2014-09-01T00:14:39.147Z
@Walter.  I also asked at http://math.stackexchange.com/questions/866836/analyze-rotation-of-satellite-orbit-due-to-transverse-acceleration.
steveOw
CC BY-SA 3.0
2014-09-01T08:16:02.010Z
Yes, there are *approximate* analytical methods (perturbation theory), valid in the limit of $K\ll1$. Perhaps you can clarify your question a bit. What's the direction of the transverse acceleration (I understand 'transverse' to mean perpendicular to the instantaneous velocity, but it's not clear whether the acceleration is in the plane of the orbit or perpendicular or a mixture).
Walter
CC BY-SA 3.0
2014-09-01T08:25:26.587Z
There is a difference between your question here and that on mathematics (and physics): here the transverse acceleration is proportional to the radial acceleration and $K$ is a dimensionless number, there the radial acceleration has no effect on the transverse acceleration and $K$ must be an acceleration (though you talk about a 'number').
Walter
CC BY-SA 3.0
2014-09-01T16:28:09.690Z
@Walter. Thanks. Apologies. This (Astronomy) version of the formula for At is what I intended. I have corrected the forumlae for At in Maths & Physics.
steveOw
CC BY-SA 3.0
2014-09-01T16:34:34.543Z
@Walter. Clarification of transverse direction: I should have emphasised that orbit is confined to a 2D plane. Then the instantaneous velocity consists of two orthogonal components (i) in the radial direction (+ away from, - towards, the Sun) and (ii) in the transverse direction (+ forwards, - backwards).  I will amend the question here.
steveOw
CC BY-SA 3.0
2014-09-01T17:17:48.560Z
What confuses me here is the word "*transverse*". I think what you mean is '*azimuthal*' (or perhaps '*tangential*'). For me, transverse means "perpendicular to the direction of motion", but that's quite different from what you meant.
Walter
CC BY-SA 3.0
2014-09-01T18:42:50.977Z
@Walter. My "transverse" direction coincides with "tangential" at perihelion & aphelion in an elliptical orbit. In a circular orbit transverse direction is the same as tangential direction. I take my lead from Calvert "Now we can find the components of the velocity. The transverse component is rdθ/dt ... The radial component is dr/dt" http://mysite.du.edu/~jcalvert/phys/orbits.htm#Coor
steveOw
CC BY-SA 3.0
2014-09-03T18:02:48.560Z
See my edited/corrected answer. The formula has changed ...
Walter
CC BY-SA 3.0
2014-09-14T11:26:36.137Z
Crossposted to http://math.stackexchange.com/q/866836/11127 and http://physics.stackexchange.com/q/123685/2451
Qmechanic
<gravity><eccentric-orbit>
Determining effect of small variable force on planetary perihelion precession
steveOw
# Determining effect of small variable force on planetary perihelion precession

(Asked by: username_0 on 2013-10-13)

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) << 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 username_1'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 username_1's Summary (in username_1'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).

**Tags:** gravity, eccentric-orbit

> **Comments:**
> **username_1**: Are you still seeking an answer?  
> **username_0**: @username_1.  Yes I am.  I have asked similar question at https://physics.stackexchange.com/questions/123685/effect-of-small-variable-transverse-force-on-satellite-orbital-precession but no solid answer received yet.  
> **username_0**: @username_1.  I also asked at http://math.stackexchange.com/questions/866836/analyze-rotation-of-satellite-orbit-due-to-transverse-acceleration.  
> **username_1**: Yes, there are *approximate* analytical methods (perturbation theory), valid in the limit of $K\ll1$. Perhaps you can clarify your question a bit. What's the direction of the transverse acceleration (I understand 'transverse' to mean perpendicular to the instantaneous velocity, but it's not clear whether the acceleration is in the plane of the orbit or perpendicular or a mixture).  
> **username_1**: There is a difference between your question here and that on mathematics (and physics): here the transverse acceleration is proportional to the radial acceleration and $K$ is a dimensionless number, there the radial acceleration has no effect on the transverse acceleration and $K$ must be an acceleration (though you talk about a 'number').  
> **username_0**: @username_1. Thanks. Apologies. This (Astronomy) version of the formula for At is what I intended. I have corrected the forumlae for At in Maths & Physics.  
> **username_0**: @username_1. Clarification of transverse direction: I should have emphasised that orbit is confined to a 2D plane. Then the instantaneous velocity consists of two orthogonal components (i) in the radial direction (+ away from, - towards, the Sun) and (ii) in the transverse direction (+ forwards, - backwards).  I will amend the question here.  
> **username_1**: What confuses me here is the word "*transverse*". I think what you mean is '*azimuthal*' (or perhaps '*tangential*'). For me, transverse means "perpendicular to the direction of motion", but that's quite different from what you meant.  
> **username_0**: @username_1. My "transverse" direction coincides with "tangential" at perihelion & aphelion in an elliptical orbit. In a circular orbit transverse direction is the same as tangential direction. I take my lead from Calvert "Now we can find the components of the velocity. The transverse component is rdθ/dt ... The radial component is dr/dt" http://mysite.du.edu/~jcalvert/phys/orbits.htm#Coor  
> **username_1**: See my edited/corrected answer. The formula has changed ...  
> **username_2**: Crossposted to http://math.stackexchange.com/q/866836/11127 and http://physics.stackexchange.com/q/123685/2451  

---

## 1 answers

### Answer by username_1 on 2014-09-01. Score: 6 | Accepted answer

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)}
$$

> **Comments:**
> **username_3**: Comments are not for extended discussion; this conversation has been [moved to chat](http://chat.stackexchange.com/rooms/16932/discussion-on-answer-by-username_1-determining-effect-of-small-variable-force-on-pla).  

---

