Here's one way to do "Minsky Motion", aka Minsky Physics. I haven't tried it with inverse square orbits, but here's a version for harmonic orbits. The Minsky method of drawing a "circle", a la Minskys & Trinskys: http://nbickford.wordpress.com/2011/04/03/the-minsky-circle-algorithm/ http://www.blurb.com/bookstore/detail/2172660?alt=Minskys+%26+Trinskys+3rd+E... x := x - delta * y ; y := y + eps * x /* _sequential_, not parallel, updates */ or equivalently, [x] [ 1 -delta ][x] [y] := [eps 1-eps*delta][y] or [x] [x] [y] := M [y], where [ 1 -delta ] M = [eps 1-eps*delta]. But we can diagonalize M with P, s.t. P^-1 M P = diag(lambda_1,lambda_2), where lambda_1, lambda_2 are the eigenvalues of M. Since det(M) = 1 = lambda_1*lambda_2, lambda_2 = 1/lambda_1. Also, trace(M) = 2-eps*delta = lambda_1+lambda_2 = lambda_1+1/lambda_1 If 0 < eps*delta < 4, then -4 < -eps*delta < 0, and -2 < 2-eps*delta < 2, and -1 < 1-eps*delta/2 < 1. So we can choose real w=acos(1-eps*delta/2), and the eigenvalues will be complex conjugates lambda_1 = cos(w)-i*sin(w) = exp(-iw), lambda_2 = cos(w)+i*sin(w) = exp(iw). However, it will be slightly more convenient to choose angle phi=w/2. Since cos(w) = 1-eps*delta/2, sin(phi) = sin(w/2) = sqrt(1-cos(w))/sqrt(2) = sqrt(1-(1-eps*delta/2))/sqrt(2) = sqrt(eps*delta/2)/sqrt(2) = sqrt(eps*delta)/2 So, phi = w/2 = asin(sqrt(eps*delta)/2), and w = 2*asin(sqrt(eps*delta)/2). Using Maxima's eigenvectors command and a lot of manipulation, [ delta delta ] P = [1-exp(-iw) 1-exp(iw)] We can now construct a real matrix [cos(-w) sin(-w)] [-sin(-w) cos(-w)] with the same eigenvalues exp(-iw), exp(iw), and [ cos(-w) sin(-w)] [exp(-iw) 0 ] Q^-1 [-sin(-w) cos(-w)] Q = [ 0 exp(iw)] [1 1] with Q = [i -i]. Combining these two efforts, we get [ 1 -delta ] [ cos(-w) sin(-w)] Q P^-1 [eps 1-delta*eps] P Q^-1 = [-sin(-w) cos(-w)] So if with compute (Q P^-1 M P Q^-1)^n, we can get the effect of n steps, i.e., [ cos(-n*w) sin(-n*w)] [-sin(-n*w) cos(-n*w)] Thus, the effect of the n-th Minsky step is M^n = P Q^-1 (Q P^-1 M P Q^-1))^n Q P^-1 In effect, we slightly skew the space, perform n iterations, and then reskew back. Slightly more perspicuously, x_n = x_0*cos(-n*w)+(x_0/2-y_0/eps)*sin(-n*w)/d and y_n = y_0*cos(-n*w)+(x_0/delta-y_0/2)*sin(-n*w)/d, where d=sqrt(1/(delta*eps)-1/4) is a real number. --------- Now instead of considering two spatial coordinates x,y, we consider instead a spatial coordinate x and its velocity -- i.e., we consider the circle in _phase space_: [x] [1 -delta][x] [v] := [eps 1-eps*delta][v] This phase plot can be interpreted as a simple harmonic oscillator consisting of a weight moving forward & back on a spring. Let k be the spring constant, and m be the mass on the end of the spring. Then m*x'' = m*v' = m*a = F = -kx, where x'' is d^2 x/dt^2. Also, conservation of energy requires that E = (1/2)*k*x^2 + (1/2)*m*v^2 = (1/2)*(k*x^2+m*v^2) If k=m in suitable units, then the phase plot of such an harmonic oscillator will be a perfect circle. Following Minsky/Trinsky, as above, we slightly skew our (digital) phase space, compute n iterations, and then reskew back to our digital phase space. Note that since the the harmonic oscillator in two dimensions is separable, we can do exactly the same thing for y, v_y, and simply include these 2 additional dimensions in our expanded phase space. For a differential equation which _isn't_ separable, we could (for n=2 objects), utilize quaternions (2 position coordinates, 2 velocity coordinates), and try to keep the _length_ of the quaternion (~ total energy) constant.