Spaces with operators notes

   Let E be a normed vector space.  Elements of L(E,E), here L(E,E) refers to a linear operator, are also called operators on E.  Let S  be a set of operators on E.  By an S-invariant subspace F we mean a subspace such that for every A in S we have AF is contained in F (i.e, if x is an element in F and A is an element in F, then Ax is an element of F). 

   An operator B is said to commute with S if AB = BA for all A an element of S.  If B commutes with S, then both the kernel of B and its image are S-invariant subspaces.  

Proof:  If x is an element of E and Bx = 0.  then ABx = BAx = 0 for all A an element of S, so the kernel of B is S-invariant.  Note: need to show if x is in ker B then Ax is also in ker B, but A(Bx) = B(Ax) = 0 means that  on the right side B(Ax) = 0 implies that Ax is in ker (B), this relation arises from the commutative relation AB = BA where the left side of our given equation provisions the A(Bx) = A(0) = 0 which ensures through associative properties that we have B(Ax) = 0 ).  We can see similarly the case for the im(B) being S-invariant.  Also recall the definition of kernel in this context meaning ker(B) = {x: Bx =0}.  It is also worth noting that the zero vector is linearly mapped to the zero vector (obvious statement I believe for normed vector spaces and linear maps...or here coupled with Banach spaces and their algebras).  

A little topology/analysis sticky for those delving into subject matter.

To prove a set is open, for any point x in a given set A defined by a given metric space, means that we need prove there exists a neighborhood N (or you might hear the terminology open ball, or open set, or other terms of a similar sort) around such point where N is contained in A...or another way of putting this, all points of N are in A, or N may be written as an open set such that N is contained in A.  

Thus in one recently read proof, the customary manner of proving this would be starting with the assumption of a given open ball of radius r, and then having shown that from the metric of any point |a-x| < r means that a is an element of A.  Here, proving that a set is open need not mean that every possible neighborhood (open set) is contained in A for a given point, but that at least one exists.  One such technique that I've seen, for example, involved to defining a radius r by the given the elements of the neighborhood metric.  Thus something like |a-x| < 1/|x^-1|  and then from this proving that a were an element of A.  In this case, set inclusion by the problem example, were implied by definition of a term of invertibility of a given element in the set, and invertibility were determined indirectly from continuity.  A presiding assumption for the given set were that |u - e| < 1 where e where a unit element (i.e. where x and e in A means that x*e = e*x = x), meant by continuity that u where in invertible (an assumption here) and thus by definition of A (in this where A were a set of invertible elements), we could see from the defined neighborhood mentioned earlier that |a*x^-1 - e| < 1 would mean in turn that the product a*x^-1 is then invertible.  Then using one other definition, namely, the product of invertible elements is in turn invertible should mean that (a*x^-1) * x  being invertible means a is invertible.

Thus a bit of multi faceted logic here working from the presumption of a given defined neighborhood which in turn leads to an element being included in a given set  which in turn proves a set being open.

  

Self note on Linear Maps and a bit of algebra

Self note on Linear maps of the following type:  Let x in R and let a be constant, and let L be a linear map, then if L(a*x/|x|) < 1, then we have a*L(x/|x|) < 1, and multiplying by |x| on both sides, we have a*|x|*L(x/|x|) < |x|.  Given linearity of L, however, we can write a*L(x/|x| +...{|x|times}...+x/|x|)= a*L(|x|*x/|x|) < |x|.  Then a*L(x) < |x|.  Any issues problems with this (normally one might expect |x| to be an integer in applying linearity here, or in other words when dealing with irrationals should there be any issues applying linearity in this manner?

Star Simulator updates

Fixed a number of issues.  One resolving properly (albeit a more minor issue) position - Center of Mass vector so that distance to the radial center of a given celestial object is computed properly in the spherical coordinate direction of position - center of mass.  Secondly, issue regarding instantaneous velocity resolved (improperly computing this from the unit vector of position - center of mass position), so that this is now computed from the properly length scale.  Thirdly, in computing big M for instantaneous velocity, instead computing for the all given celestial masses whose radial magnitude is such that r_s > r where r is the distance of any star less than a given star s.

Updates finished.

At initialization portion of code, you can change the distribution of mass along the z axis by changing a numeric divisor where approaching 1 means that a given distribution is expressed by the original z axis parametric random distribution given.  You might notice, that changing this value, with greater mass distribution of stars outside of  a quasi 2 dimensional rotational plane, and having added stars into a given system, tends to lead to a non zero momentum  barycentric system.  While those systems tending towards flattened rotational galactic planes to be more stable.  As opposed to a system that has no initial linear momentum, also it would be noticed that the frequency of stars gobbling other stars is less so.  

Rotational plane star simulation

Okay I from a past class in Astronomy.  I grabbed this program and tweaked it a bit.  I managed to add initial momentum into a given system, and I were attempting to simulate something of a galactic rotational plane (without expansion).    Several variables worth noting for changing the simulator.  L changes the distance between stars in general, while expanding the field of view.  Its worth noting that you'd likely in order to see the stars better when say expanding out by or in by powers of ten correspondingly changing the
variable pertaining the visible radius set for stars.

This is found on the line

Stars = Stars+[sphere(pos=(x,y,z), radius=10*r, color=col0)]

While I've tried using something of a centripetal acceleration for tangential motion of stars related to a given center of mass, the effort is off or even in this case, setting a larger scale mass (with zero initial momentum),
relative to a given star mass, I've still encountered some issues, namely, an expansion in the overall system.  I've managed to tweak, however, initial momentum of stars in general so that rotation is about as stable as I've seen.  My system is a bit slow in processing, so I couldn't remotely get close to something like our own galaxy in the number of stars, and likely one would be dealing with much greater length scales here, relative to the given length scales presented in L.
Adjusting the parameter cmass/6e10: vmags.append([vmag(cmass/6e10,rmag1)])

from visual import *
from time import clock
from random import random,uniform
import numpy
import math

# Sets up a bunch of masses and lets them move under the influence of their
# mutual gravity.
# Masses are put in random positions within a sphere originally. An additional
# force is applied which simulates an equal density distributed over the
# rest of the universe outside this initial sphere.


def randomdirection():
      # Generates random direction on sky.
      # RA (in radians)
      # Dec (in radians)
      ra = 2.0*math.pi*random()
      dec = math.acos(2.0*random()-1.0)-0.5*math.pi
      return ra,dec

def ranvec():
      # Generates a randomly orientated unit vector.
      theta, phi = randomdirection()
      z = math.sin(phi)
      x = math.cos(phi)*math.sin(theta)
      y = math.cos(phi)*math.cos(theta)
      vec = numpy.array([x,y,z])
      return vec


# Stars interacting gravitationally
# Program uses numpy arrays for high speed computations


Nstars = 100  # change this to have more or fewer stars

G = 6.7e-11 # Universal gravitational constant

# Typical values
Msun = 2E30
Rsun = 6.96e8
Rtrail = 2e8
L = 4e11
vsun = 0.9*sqrt(G*Msun/Rsun)
h0 = 2.5e-6 # Hubble's constant - expansion rate 2.5e-6 is good.

scene = display(title="Stars", width=1300, height=740,
                range=L, forward=(-1,-1,-1))


Stars = []
poslist = []
plist = []
mlist = []
rlist = []
p0 = 0.0*Msun*100000.0

for i in range(Nstars):
    vec = L*ranvec()*(random()**0.3333)
    x = vec[0]
    y = vec[1]
    z = vec[2]/1e4
    r = Rsun
    col0 = (uniform(0.7,1.0),uniform(0.7,1.0),
                  uniform(0.7,1.0))
    Stars = Stars+[sphere(pos=(x,y,z), radius=10*r, color=col0)]
    mass = Msun
    px = p0*uniform(-1,1)
    py = p0*uniform(-1,1)
    pz = p0*uniform(-1,1)
    poslist.append((x,y,z))
    plist.append((px,py,pz))
    mlist.append(mass)
    rlist.append(r)

pos = array(poslist)
p = array(plist)
m = array(mlist)
m.shape = (Nstars,1) # Numeric Python: (1 by Nstars) vs. (Nstars by 1)
radius = array(rlist)

vcm = sum(p)/sum(m) # velocity of center of mass
def centerofmass(m, pos):
    # find center of mass
##    print(m.shape)
##    print(pos.shape)
##    print(pos*m)
##    
##    print(sum(pos*m,axis=0))
    R = sum(pos*m,axis=0)/sum(m)
    return R

def rfromRvec(pos,R):
    return pos-R

def unitnorm(r2):
    v = []
    for x in r2:
        v.append([numpy.linalg.norm(x)])
    v = array(v)
    return r2/v

def tangentfromr2(r2):
      def orthogonalangle(x,y):
            return math.pi/2 + math.atan(y/x)
      def findtancoords(rxpyp,oalpha):
            return [rxpyp*math.cos(oalpha),
                    rxpyp*math.sin(oalpha)]
      def findphi(rxy,z):
            return math.atan(z/rxy)
      def findrxy(x,y):
            return (x**2+y**2)**.5
      def findrxyz(x,y,z):
            return (x**2+y**2+z**2)**.5
      def rotationphi(phi):
            m = [(math.cos(phi),math.sin(phi)),
                 (-math.sin(phi),math.cos(phi))]
            return numpy.array(m)
      def revrotationphi(phi):
            m = [(math.cos(phi),-math.sin(phi)),
                 (math.sin(phi),math.cos(phi))]
            return numpy.array(m)      
      def rotatecoords(x,z,R):
            xzcrd = [x,z]
            XZ = numpy.array(xzcrd)
            rXZ = numpy.dot(R,XZ)
            return rXZ

      
# to find tangent one solution, involves using rotated coordinate system
# on xz plane by phi where phi is represented by the angle of r2.
#  Here we then have under rotated coordinates a tangent vector given
# by a perpendicular on the 2 dimensional plane by pi/2, then covert
# computed unit tangent vector back to our original coordinate system.
# a little easier despite the fuss of coordinate transformations since
# computing a tangent on a 2 dimensional plane is a little less cumbersome.

      #first find r for xy plane
      x=r2[0];y=r2[1];z=r2[2]
      rxy = findrxy(x,y)
      theta = findphi(x,y)
      if theta < -math.pi:
            print(theta)

      phi = findphi(rxy,z)
      if phi <0:
            print(phi)
      Rxy = rotationphi(theta) 
      Rxz = rotationphi(phi)
      revRxz = revrotationphi(phi)
      revRxy = revrotationphi(theta)
      #phi positive is 1rst or 3rd quad
      #theta positive is 1rst or 3rd quad
      XpYp = rotatecoords(x,y,Rxy)
      #z = zp
      xp = XpYp[0]; yp = XpYp[1]; zp = z;
      XppZpp = rotatecoords(xp,zp,Rxz)
      xpp = XppZpp[0]; zpp = XppZpp[1]; ypp = yp
      # under coordinate rotations ypp = yp,
      
##      orthetap = 1*math.pi/2
      #here a twice rotated coordinate system puts the zpp at zero, and likewise ypp should be zero)
      rxppypp = findrxy(xpp,ypp)
      typp, txpp = (xpp,ypp)
##      tangentcoords = findtancoords(rxppypp, orthetap)
##      txpp,typp = tangentcoords
      if typp < 0:
            print(typp)
##      typp = abs(typp); txpp = abs(txpp)
      #convert txp, zp back to original coordinates, typ = ty
      Txpzp = rotatecoords(txpp,zpp,revRxz)
      txp = Txpzp[0]; typ = typp; tzp = Txpzp[1]
      Txy = rotatecoords(txp,typ,revRxy)
      tx = Txy[0]; ty = Txy[1]; tz=tzp
      return [tx,ty,tz]

def vmag(M,rmag):
      # M is the total mass, and rmag is the distance from center of total mass M to star
      G = 6.7e-11
      return (G*M/rmag)**.5

def angles(pos):
      #find spherical coordinate angles
      c = pos[:,1]/pos[:,0]
      d = (pos[:,0]**2 + pos[:,1]**2)**.5
      e = pos[:,2]/d
      f = []; g = [];
      i = 0
      for x in c:
            f.append([math.atan(x),math.atan(e[i])])
            i += 1

      return numpy.array(f)

def finddcoords(rxy,oalpha):
      r2 = []
      for oalphai in oalpha:
            theta,phi = oalphai
            rx = rxy*math.sin(theta)*math.cos(phi)
            ry = rxy*math.sin(theta)*math.sin(phi)
            rz = rxy*math.cos(theta)
            r2.append([rx,ry,rz])
      return numpy.array(r2)

def getCMlist(r2,m):
      magnr2 = (r2[:,0]**2+r2[:,1]**2+r2[:,2]**2)**.5
      i = 0;Mlist =[];
      for m2 in magnr2:
            n=0;M=0
            for m2a in magnr2:
                  if not i == n:
                        if m2a < m2 or m2a == m2:
                              M += m[n]
                  n+=1
            Mlist.append(M)
            i+=1
      return numpy.array(Mlist)
                                    
      
m[0] = sum(m)/1
cmass = sum(m)
R = centerofmass(m, pos)
print('center of mass:')
print(R)
pos[0] = R
R = centerofmass(m, pos)
M = sum(m)
#Rsun in the direction of r
ang = angles(pos)
rSunarr = finddcoords(Rsun,ang)
r2 = rfromRvec(pos,R) + rSunarr
Mlist = getCMlist(r2,m)
print(r2)
r2n = unitnorm(r2)
print(r2n)
tpos = []
vmags = []
itr = 0
print('mlist:')
print(Mlist)
for r2i in r2n:
      txyz = tangentfromr2(r2i)
      tpos.append(txyz)
      x = r2[itr,0]; y=r2[itr,1]; z=r2[itr,2];
      rmag1 = (x**2+y**2+z**2)**.5
##      vmags.append([vmag(cmass/6e10,rmag1)])
      vmags.append([vmag(Mlist[itr],rmag1)])
      itr+=1
tposa = array(tpos)
print(tposa)
vmagsa = array(vmags)
print(vmags)
vs = tposa*vmagsa
print(vs)

p = p-m*vcm # make total initial momentum equal zero
p = vs*Msun
p[0] = [0,0,0]



print(p)
dt = 100.0
Nsteps = 0
pos = pos+(p/m)*(dt/2.) # initial half-step
time = clock()
Nhits = 0

while 1:
    rate(10000)

    L *= 1.0+h0*dt
    con = 0.9*G*Nstars*Msun/(L*L*L)# strength of force to allow for external mass
    
    # Compute all forces on all stars
    try:  # numpy
        r = pos-pos[:,newaxis] # all pairs of star-to-star vectors
        for n in range(Nstars):
            r[n,n] = 1e6  # otherwise the self-forces are infinite
        rmag = sqrt(add.reduce(r*r,-1)) # star-to-star scalar distances
        hit = less_equal(rmag,radius+radius[:,newaxis])-identity(Nstars)
        hitlist = sort(nonzero(hit.flat)[0]).tolist() # 1,2 encoded as 1*Nstars+2
        F = G*m*m[:,newaxis]*r/rmag[:,:,newaxis]**3 # all force pairs
    except: # old Numeric
        r = pos-pos[:,newaxis] # all pairs of star-to-star vectors
        for n in range(Nstars):
            r[n,n] = 1e6  # otherwise the self-forces are infinite
        rmag = sqrt(add.reduce(r*r,-1)) # star-to-star scalar distances
        hit = less_equal(rmag,radius+radius[:,newaxis])-identity(Nstars)
        hitlist = sort(nonzero(hit.flat)) # 1,2 encoded as 1*Nstars+2
        F = G*m*m[:,newaxis]*r/rmag[:,:,newaxis]**3 # all force pairs
        
    for n in range(Nstars):
        F[n,n] = 0  # no self-forces
    p = p+sum(F,1)*dt+pos*con*dt*m

    # Having updated all momenta, now update all positions         
    pos = pos+(p/m)*dt

    # Expand universe
    #pos *= 1.0+h0*dt

    # Update positions of display objects; add trail
    for i in range(Nstars):
        Stars[i].pos = pos[i]

    # If any collisions took place, merge those stars
    for ij in hitlist:
        i, j = divmod(ij,Nstars) # decode star pair
        if not Stars[i].visible: continue
        if not Stars[j].visible: continue
        # m[i] is a one-element list, e.g. [6e30]
        # m[i,0] is an ordinary number, e.g. 6e30
        newpos = (pos[i]*m[i,0]+pos[j]*m[j,0])/(m[i,0]+m[j,0])
        newmass = m[i,0]+m[j,0]
        newp = p[i]+p[j]
        newradius = Rsun*((newmass/Msun)**(1./3.))
        iset, jset = i, j
        if radius[j] > radius[i]:
            iset, jset = j, i
        Stars[iset].radius = newradius
        m[iset,0] = newmass
        pos[iset] = newpos
        p[iset] = newp
        Stars[jset].visible = 0
        p[jset] = vector(0,0,0)
        m[jset,0] = Msun*1E-30  # give it a tiny mass
        Nhits = Nhits+1
        pos[jset] = (10.*L*Nhits, 0, 0) # put it far away

    Nsteps += 1

Specific solution method for computing a tangent vector in a 3 dimensional case for a given plane of rotation

Of course, it seems dot product formulation u dot v = 0 seems straight forward enough in computing an orthogonal vector.  The problem is that there can be a lot of different orientations for a given orthogonal vector, and thus many different solutions to the same problem (I'd have an infinite number of solutions, for instance, rotated around the plane of r and z defined by angle phi on the set of reals with 2 solutions matching the given plane).  So one approach that I've considered, not sure if its a bit lengthy but at least it provides a bit of user control in the direction choosing aspect of an orthogonal vector.  Is in using a bit of trigonometry.  Lets say, in a given coordinate system we knew the direction of r(x,y) by its angle theta, and on the z plane its angle between r and z given by phi.  Then we could do a bit of work in three dimensions with these angles by adding say 90 degrees to a given angle and then recomputing for the new given direction of r or z, a given set (x,y,z) coordinates.  Another solution, that I have found is using a bit of linear algebra and using a rotation of coordinates system, where rotating our coordinate plane xz to match the angle of phi.  Where under the new coordinate system we have the vector now occupying a two dimensional plane.  In this case finding an orthogonal vector, at least in theory, sounds a bit easier, here we have our coordinate plane matched to the angle between r and z which means that we aren't merely computing a tangent that happenstance is incidentally contacting the surface of such plane of rotation, but is matched to the plane that we have in mind.  In the case of a two dimension tangent computation, we'd likely have 2 possible solutions for such plane, and thus we have a bit of measured control over potential solutions to be yielded here.   Once computing the new angle of  r(x'',y'') I'd call theta' in the rotated coordinate axis, we can reconvert the tangent coordinates from the rotated coordinate system back into our original coordinate system, and voila we have a tangent vector that is situated to a given rotational plane.  Here's a bit of python code expressing this trigonometric solution.  You'd need to have numpy and python 2.6 installed on your system for this

def unitnorm(r2):
    v = []
    for x in r2:
        v.append([numpy.linalg.norm(x)])
    v = array(v)
    return r2/v

def tangentfromr2(r2):
      def orthogonalangle(x,y):
            return math.pi/2 + math.atan(y/x)
      def findtancoords(rxpyp,oalpha):
            return [rxpyp*math.cos(oalpha),
                    rxpyp*math.sin(oalpha)]
      def findphi(rxy,z):
            return math.atan(z/rxy)
      def findrxy(x,y):
            return (x**2+y**2)**.5
      def findrxyz(x,y,z):
            return (x**2+y**2+z**2)**.5
      def rotationphi(phi):
            m = [(math.cos(phi),math.sin(phi)),
                 (-math.sin(phi),math.cos(phi))]
            return numpy.array(m)
      def revrotationphi(phi):
            m = [(math.cos(phi),-math.sin(phi)),
                 (math.sin(phi),math.cos(phi))]
            return numpy.array(m)      
      def rotatecoords(x,z,R):
            xzcrd = [x,z]
            XZ = numpy.array(xzcrd)
            rXZ = numpy.dot(R,XZ)
            return rXZ

So basically under a twice rotated coordinate system, we match for the angle theta on the xy plane to that of r, and then rotate to the x'z' plane to match the angle phi for coordinates x'',y'', z'' . This in turn should lead us to a vector whose coordinate is (x'', 0, 0), we can then simple compute on a given plane the orthogonal by computing the new direction of r at 90 degrees which should leave a component (0,y'',0) and then twice rotate back into our original coordinate system.

If you wanted to change the respective rotational plane you could rotate the coordinate plane from here on the y''z'' plane in any given orientation and then use the same process (90, or -90 degree) computation for the new vector r in the x''', y''', z''' coordinate system, and then reverse the process for coordinate transformations back into x, y, z.  Thus this method provides a solution method for defining our given rotational plane and computing tangents on such a plane.

Hmm..Saturnine...