Numerical error in implementation of iterative algorithm

I am trying to implement (in Python for now) low thrust orbit propagation for spacecraft using universal variables. For a given central body with the gravitational parameter $\mu$ and an orbit with the semi-major axis $a$ and the initial position $\vec{r}_0$ and velocity $\vec{v}_0$ at $t = t_0$ the position for a given time is given by: $$\vec{r} = \vec{r}_0 \ f(s) + \vec{v}_0 \ g(s)$$ and the velocity by $$\vec{v} = \vec{r}_0\ \dot{f}(s) + \vec{v}_0\ \dot{g}(s)$$

Where $$f(s) = 1-\left(\frac{\mu}{|\vec{r}_0|}\right) s^2 c_2(\alpha s^2)$$ $$g(s) = t-t_0-\mu s^3c_3(\alpha s^2)$$

With $\alpha=\frac{\mu}{a}$, $c_n$ the n-th Stumpff function and $s$ the solution of $$t-t_0=|\vec{r}_0| s \ c_1(\alpha s^2) + |\vec{r}_0| s^2 \ c_2(\alpha s^2) + \mu s^3 \ c_3(\alpha s^3)$$

I am calculating $s$ by using Newton's Method and everything works fine for scenarios without thrust. The orbit is elliptical and closed, parabolic for higher initial velocities, if I choose a $\vec{r}_0$, $\vec{v}_0$ and $t_0$ and propagate $t$ forward.

For cases with thrust where I would change the velocity in every iteration

$$\vec{v}_{k} = \vec{v}_{k-1} + \vec{a}_{k} \cdot \Delta t$$

it is necessary to update $\vec{r}_0$ and $\vec{v}_0$ every time, propagate $\Delta t$ forward and repeat instead of choosing initial conditions and just propagating $t$ forward. In preparation for that I chose not to add any velocity so both approaches should yield the same results.

Here is my Problem: Even after the first iteration the differences are large and the error increases fast blue: fixed initial conditions, propagating $t$

red: propagating initial conditions, fixed $\Delta t$

I thought I could solve this by using a more advanced integration algorithm like Runge-Kutta and tried to calculate every iteration in as few steps as possible but I was not able to transform the equations so that I could use Runge-Kutta (since it is not an ODE) and reducing the steps did not help at all.

Can someone help to fix this or give any hints to why this error is so large? Thank you in advance!

P.S.: I already ask this question on the physics stack exchange and they advised me to ask it here.

EDIT: Adding the code: Pastebin

• If decreasing $\Delta t$ or increasing your precision for $s$ does not resolve the matter then it's presumably an implementation issue in which case you'll need to post your code for help. – lemon Aug 18 '16 at 5:47
• I have added the code – Ganymed_ Aug 18 '16 at 17:11
• And yes, of course decreasing $\Delta t$ helps but still energy is not conserved and the orbit changes overtime – Ganymed_ Aug 18 '16 at 17:26
• While I can't spot the error with a cursory look at the code, it seems clear to me from your results that the error is with your evaluation of the new velocity... – lemon Aug 18 '16 at 19:18