Determining effect of small variable force on planetary perihelion precession
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)}
$$