While the equations
v(n+1) = v(n) + a
d(n+1) = d(n) + (v(n+1) + v(n)) / 2
- d(n) is position at time n;
- v(n) is velocity at time n; and
- a is acceleration in distance units per frame per frame
are arithmetically (and physically) correct, they are computationally problematic.
The faster your frame rate becomes the larger the computational error arising from floating point approximation becomes as a fraction of your values for velocity, acceleration and time. Then the further out you extrapolate these approximations from the last change in acceleration, the greater the accumulated error in position and velocity becomes.
The problems arise in attempting to calculate intersection of targeting vectors, or predict future positions.
A better computation is to directly calculate position and velocity from the initial values at the last acceleration change and the elapsed time as:
v(n) = v(0) + a * n
d(n) = d(0) + (a * n * n) / 2
where v(0) and d(0) are respectively the position and velocity at the time of the last acceleration change, and n is the number of frame time-units since the last acceleration change.
To compute the gravitational force between two objects one uses Newton's gravitational law:
F = m1 * m2 * G / (r * r)
and the gravitational acceleration experienced by each mass due to the other will be
g = M * G / (r * r)
where G is the Gravitational Constant
Note that the gravitational acceleration experienced by each mass is only proportional to the mass of the other object, not it's own.
Generalizing the homily above on computational accuracy and precision, with non-constant but calculable acceleration such as in this case, it would be better to use the appropriate Calculus formulae to calculate position and velocity directly at each frame than to calculate each indirectly from the immediately previous values.