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+…
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.
