I should have pointed out that the Trinsky analysis of Minsky essentially produces Pythagorean triples/Gaussian primes for certain rational numbers; e.g., if d=1, delta*eps = 4/(1+4*d^2) = 4/(1+4*1^2)=4/5, and w = 2*asin(sqrt(4/5)/2) = 2*asin(1/sqrt(5)) = asin(4/5) -- i.e., we get a 3,4,5 triangle. When we compute arbitrarily high powers ((3+4*i)/5)^n, we conserve the length (i.e., energy) of the system, but we don't actually get an "orbit", since Q(i) is a field of characteristic 0. I now need to analyze the much harder case where we execute the _floor(x)_ and/or _round(x)_ function at each Minsky step. This system is the one actually discovered by Minsky, and produces a slightly blurry circle/ellipse, because it touches on a number of close by integer points and may require several "cycles" before it literally repeats.

From: Henry Baker <hbaker1@pipeline.com>

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.