Sunday, May 29, 2016

Space Elevator Calculations

In this post we use the equations derived in the previous posts and allow the user to plug in some numbers.
In particular, when a cable thickness has been chosen to have the optimal thickness, we determine the mass of satellite required to achieve the given tension at the earth's surface.
We also work out the total mass of the cable.

Inputs
Cable Specific Strength M Yuri
Safety factor -
Tension at earth's surface Newtons
Satellite height above earth km

Outputs
Total Cable mass - kg
Satellite Mass - kg
Tension at Satellite - Newtons

Saturday, May 7, 2016

Lifting a cable above the equator

Suppose you were at the equator and decided to hold a cable at some height above the earth's surface. The cable dangles vertically down below, touching the earth, but it is not anchored there. As you go higher and higher, you would stay above the same point on the earth. I.e. the top of the cable would be geostationary as the earth spins.
If the cable were not made of a very strong material then as you tried to go higher, it would eventually break.
It turns out that if you have a cable of specific strength 48.5 M Yuri, then you'd be able to reach all the way up to a geo-stationary satellite (35,797km up).
If you reach that height, then you'll be able to go on out to about 143,937 km without requiring a stronger cable.
Something a tad surprising ( to me at least) happens 143,937 km up. You will no longer need to hold the top end. The tension at both the top and the bottom will be zero. The maximum tension will be at the height of a geostationary satellite orbit.
A cable that long, will be able to support itself. The gravitational attraction is exactly the required centrifical force to maintain a circular orbit.
If you then attempted to use a longer cable and go higher, the cable would lift off the ground. In other words the net force would be upwards. So if you want to go out further, then you'd need to anchor one end to the surface of the earth.

I've made lots of statement there, now lets start to justify some of that mathematically.
This post is a follow on from the previous one: Polar Escape Strength Threshold .
We'll use the same variables as defined there.
The cable is vertical and as the earth spins, the top of the cable remains over the same point on the equator. So the cable is geostationary and hence the net force acting on a small section of cable is the centrifical force: \[ \delta m\ \omega_s^2 r = \frac{G M_e \delta m}{r^2} - \delta T \] but we know: \[\delta m = \lambda \ \delta r\] and so: \[ \delta T = \delta r \left( \frac{G M_e \lambda}{r^2} - \lambda \omega_s^2 r\right) \] We can integrate that and use the fact that the tension in the cable at the earth's surface is zero. \(T(r_e) = 0 \).
We find:
\[ T(R) = G M_e \lambda \left( \frac{1}{r_e} - \frac{1}{R} \right) + \frac{\lambda \omega_s^2}{2} \left( r_e^2 - R^2 \right) \ \ \ \ \ \ \ \ \ (eqn 1) \] We define the specific strength Y to be the tension at breaking point divided by the mass per unit length: \[ Y = \frac{T(R_b)}{\lambda} = G M_e \left( \frac{1}{r_e} - \frac{1}{R_b} \right) + \frac{ \omega_s^2}{2} \left( r_e^2 - R_b^2 \right) \]
We now return to eqn 1 above. We can differentiate that with respect to R to find where is the maximum tension: \[ 0=\frac{\partial T}{\partial R}= \frac{G M_e \ \lambda}{R^2}- \lambda \ \omega_s^2 R \] We find that the maximum tension occurs at the geostationary orbit: \[ R_g = \left( \frac{G M_e}{\omega_s^2} \right)^{1/3} \] The threshold specific strength that would enable the cable to reach this max tension is: \[ Y_E = GM_e \left( \frac{1}{r_e} - \frac{1}{R_g}\right) + \frac{\omega_s^2}{2} \left( r_e^2 - R_g^2\right) \approx 48.5 M Yuri \] Returning to eqn 1, we see that \(T(R)=0\) when \(R=r_e\), in other words the tension is zero at the earth's surface. The tension reaches a maximum at the orbit of a geostationary satellite. It then decreases again. But we can find R such that \(T(R)=0\) \[ 0 = G M_e \left( \frac{1}{r_e} - \frac{1}{R_s} \right) + \frac{ \omega_s^2}{2} \left( r_e^2 - R_s^2 \right) \] That implies: \[ 0 = R_s^3 \ \frac{ \omega_s^2} {2} - R_s \left( \frac{ \omega_s^2 \ r_e^2 }{2} + \frac{G M_e }{r_e} \right) + G M_e \] That is cubic in \(R_s\). We solve that numerically and we find that the distance from the centre of the earth to the point on the cable where the tension is zero is: \[R_s \approx 1.5030 \times 10^8 m\] which 143,937 km from the surface of the earth.
A cable which is that long and is placed vertically above a point on the equator has zero tension at both ends. It neither drifts off into space nor falls to the earth. It just remains where it is.

Returning once again to my (current) favorite equation i.e. eqn 1, we see that it implies that for \(R>R_s\) the tension becomes negative. But we can't allow negative tension since we can't push with our cable. Without being held down the cable would rise up. So we would need an anchor at the surface of the earth. With that anchor we no long would have the conditon \(T(r_e) =0 \) . That was used in the derivation of eqn 1. And so eqn 1 would not be valid.
What we find is that as the cable gets longer beyond 143,937 km, the tension would keep increasing. An infintely long cable would require infinite tension at the anchor.

Polar Escape Strength Threshold

Suppose we had a cable with mass per unit length \(\lambda\).
So the mass of a small section of length \(l\) is: \[m \ = \lambda \ l \] If the cable is hanging vertically in a constant gravitational field \(g\), then the tension at the top will be the \[T = m \ g = \lambda \ l \ g \] We define the specific strength to be the tension at breaking point divided by the mass per unit length: \[Y=\frac{T_b }{\lambda} \] This is measured in units of: \[[Yuri] = \frac{[Newton] [ meter]} {[kilogram]^2} \] With constant \(g\), a cable with specific strength \(Y\) will be able to support up to its breaking length: \[l_b =\frac{T_b }{\lambda \ g} = \frac{Y}{g}\]
Now how about the case when the gravitational follows a one over r squared law. Where r is the distance from the centre of the earth. Lets start by considering a vertical cable that isn't moving: A small section from \(r\) to \((r + \delta r)\) has forces acting up and down. If it is not accelerating, then the forces are balanced: \[F_{up}=F_{down}\] \[T(r + \ \delta r )= T(r)+ \frac{G M_e \delta m }{r^2} \] But we know the mass: \[\delta m = \lambda \ \delta r\] So \[\delta T = \frac{G M_e \ \lambda \ \delta r }{r^2} \] We can integrate both sides of that for a cable that goes from the earth's surface \(r=r_e\) to \(r=R=r_e + h \) a height h above the earth.
We note that the end that touches the surface is not anchored. So the tension: \[T(r_e)=0\] So \[T(R)=G M_e \ \lambda \left( \frac{1}{r_e} - \frac{1}{R} \right)\] Suppose the cable breaks at \(R_b\), returning to the definition of the specific strength: \[Y=\frac{T(R_b)}{\lambda} \] Thus: \[Y=G M_e \ \left( \frac{1}{r_e} - \frac{1}{R_b} \right)\] rearranging we find: \[R_b = \left( \frac{1}{r_e} - \frac{Y}{G M_e}\right)^{-1}\]
So the breaking height above the earth's surface: \[h_b = R_b - r_e = \left( \frac{1}{r_e} - \frac{Y}{G M_e}\right)^{-1} - r_e\] Now when: \[\frac{1}{r_e} - \frac{Y}{G M_e} = 0 \] we can reach an infinite height! In other words the cable can completely escape the earth's gravitational field.
The polar escape strength threshold is: \[Y_p = \frac{G M_e}{r_e} \approx 62.6 M Yuri \] We use the word polar here since is most relevant for a cable above either the north or south pole. If we were to have a geostationary cable above the pole, then we could ignore the spin of the earth. On the other hand for a geostationary cable above the equator, the centrifical acceleration should not be ignored. I'll show the details of how to deal with that in my next post.

Monday, May 2, 2016

Sky Hooks

If we want to make it easy(ish) to get into space, then it would be very helpful to have a sky hook. Once we have that in place we can build an elevator.
Lets look into the physics / maths of how it could be done.
Suppose a small object of mass m is travelling in a cirle round a sphere. The object is a distance R from the centre of the sphere.
In vector notation, write its position as: \[\textbf{p}=R\left( \begin{array}{c} \sin(\omega t) \\ \cos(\omega t) \end{array} \right)\] We can then differentiate that twice with respect to time and we find the acceleration: \[\textbf{a}=\frac{\partial^2\textbf{p}}{\partial t^2}=-\omega^2 R \left( \begin{array}{c} \sin(\omega t) \\ \cos(\omega t) \end{array} \right) = - \omega^2 \textbf{R} \]
So, when an object is moving in a circle, there is an acceleration towards the centre of that circle.
But we know that the force: \[\textbf{F}=m \textbf{a}\] So if we want to maintain the object rotating about the centre, then we need apply a force toward the centre of magnitude: \[F=m \omega^2 R\] That is called the centrifical force.
For a satellite travelling round a planet, there is a gravitational attraction: \[F_g=\frac{G M_e m}{R^2} \] Now, if the gravitational attraction is equal to the centrifical force, then the satellite can travel in a circular orbit: \[ m \omega^2 R = \frac{G M_e m}{R^2} \] We can solve for \(R\): \[R = \left( \frac{G M_e}{\omega^2} \right)^\frac{1}{3}\]
For a geostationary orbit, as the earth spins on its axis, the satellite will remain above the same point on the equator, we set: \[ \omega_s \ n_d s_d =2 \pi n_r \] where the number of seconds in a day: \(s_d = 24 * 3600\)
the number of days in a year \(n_d\approx365.2422\)
the number of revolutions of the earth about its own axis in a year \(n_r=n_d+1\)
(The +1 there comes from the fact that the earth is also revolving round the sun, for details see the Sidereal time article in Wikipedia).
And we find: \[\omega_s \approx 7.29211 \times 10^{-5} \ rad . s^{-1}\] and plugging that into the equation for \(R\) we find: \[R_g \approx 4.2164 \times 10^7 m = 42,164 km \] That's the distance from the centre of the earth to a geostationary satellite.
It is 35,786 km above sea level.

Suppose now we were to consider a different geostationary satellite, this time with a cable attached. The other end of the cable is anchored to the ground.
Now suppose the cable is perfectly vertical and there is a tension \(T(r)\) in the cable at r. The net force on the satellite should be the centrifical force: \[T (R) + \frac{G M_e M_s}{R^2}= M_s \omega_s^2 R \ \ \ \ \ \ \ \ \ \ \ \ (eqn 1)\] Now we consider a small section from \(r\) to \(r+\delta r\) above the centre of the earth.
Suppose the cable had constant mass per unit length \(\lambda\). The mass of the small section is: \[\delta m = \lambda \delta r\] The gravitational force on that section is: \[\frac{G M_e \lambda \delta r}{r^2}\] We write the net force on the small section from \(r\) to \(r+\delta r\) due to the change in tension as: \(\delta T(r)\).
For the cable section to remain in circular orbit we need the net force to be the centrifical force. Hence: \[\frac{G M_e \lambda \delta r}{r^2} - \delta T(r) = \lambda \delta r \omega_s^2 r\] \[\delta T(r) = \frac{G M_e \lambda \delta r }{r^2} - \lambda \delta r\omega_s^2 r\] \[ T(R) - T(r) = \int^R_r dT(r') = G M_e \lambda \left( \frac{1}{r} - \frac{1}{R} \right) + \frac{1}{2} \lambda \omega_s^2 \left( r^2 - R^2 \right) \] But we already have an expression for \(T(R)\) above (eqn 1). So, after a bit of re-arranging we find the tension in the cable to be: \[ T(r) = G M_e \left( \frac{\lambda}{R} - \frac{\lambda}{r} -\frac{M_s}{R^2} \right) + \frac{1}{2} \lambda \omega_s^2 \left( R^2 - r^2 \right) +M_s \omega_s^2 R \ \ \ \ \ ( eqn 2)\] Let \(r_m\) be the position of the maximum tension in the cable. \[0=\frac{\partial T}{\partial r}(r_m)\] We find that occurs at the orbit of an untethered geostationary satellite: \[r_m = R_g = \left( \frac{G M_e}{\omega_s^2} \right)^\frac{1}{3}\] You can't push with our rope!
The cable has tension but does not offer vertical support.
Mathematically we write that as: \[T(r) \geq 0\ \ \ \ \ \ \ \ \forall\ r : r_e \leq r \leq R \] To support the mass of the cable, the satellite must be further out than \(R_g\)
We can use eqn 2 to determine the tension in the cable at the earth's surface \( T(r_e)\).
That's the maximum force the climber is allowed exert so as not to pull down the satellite.

Thin in the right places

Suppose the breaking tension of a cable is \(T_b\) then we define the specific strength \(Y\) to be: \[Y = \frac{T_b}{\lambda} \] We define an optimal non-uniform cable to be one where at every point it is just thick enough to safely withstand the tension at that point.
Mathematically we write that as: \[ \lambda(T ) = \frac{T\ \alpha}{Y} \ \ \ \ \ \ \ \ \ (eqn 3) \] where our \(\alpha\) is the safety factor with: \[1 \leq \alpha \] With a little algebra we find: \[\alpha = \frac{T_b}{T}\] and so we know that \[T \leq T_b\] in other words the tension in the cable is at or below the tension at which the cable breaks. A higher \(\alpha\) is safer.
But now the mass of a small section from \(r \rightarrow r + \delta r\) is: \[\delta m = \delta r \lambda \left( T (r) \right) = \delta r \ \frac{\alpha T(r)}{Y}\] Just as before, to keep the cable section in a circular orbit the net force should equal the centrifical force. The net force is made up of a combination of the change in tension and the gravitational force. \[T (r) - T (r + \delta r) + \frac{G M_e \ \delta r \ \alpha T (r)}{Y \ r^2} = \frac{\omega_s^2 r \ \delta r \ \alpha T (r)}{Y}\] We can rearrange that: \[ \frac{\delta T}{T} =\delta r \left[ \frac{G M_e \alpha}{Y r^2} - \frac{\omega_s^2 r \alpha}{Y} \right] \] We now integrate both sides between \(r_1\) and \(r_2\) where \[ r_e \leq r_1 < r_2 \leq R \] And we find: \[ ln \left( T(r_2) \right) -ln \left( T(r_1) \right) =G M_e \frac{\alpha}{Y} \left[ \frac{1}{r_1} - \frac{1}{r_2}\right] + \frac{\omega_s^2 \alpha}{2 Y} \left( r_1^2 - r_2^2 \right) \] and we can write that as: \[ \frac{T(r_2)}{T(r_1)} =exp \left( G M_e \frac{\alpha}{Y} \left[ \frac{1}{r_1} - \frac{1}{r_2}\right] + \frac{\omega_s^2 \alpha}{2 Y} \left( r_1^2 - r_2^2 \right) \right) \ \ \ \ \ (eqn 4) \] Now bear in mind that when deriving that equation we have used an optimally thin cable, (see equation 3). At all points the cable mass per unit length has been chosen to be just thick enough to be safe.
Now if we know the tension at any one spot from ground all the way up to the satellite, then we can use equation 4 to work at the tension at all other locations along the cable. We could specify the tension at the ground anchor and then work out the tension all the way up to the satellite.
We consider 3 variables:
\(M_s\): the mass of the satellite.
\(R\): the distance of the satellite from the centre of the earth.
\(T(r_e)\): the tension at the anchor on the surface of the earth.

Now if we have any two of those we can work out the third.
For example we could specify the distance of the satellite and the tension at the ground anchor and we can work out the mass of the satellite required.
To do that, we would use equation 4 to work out the tension in the cable when it meets the satellite: \(T(R)\).
But we also know that at the satellite, the centrifical force is equal to the gravity plus the tension, so: \[ M_s \omega_s^2 R = T(R) + \frac{G M_e M_s}{R^2} \ \ \ \ (eqn 5)\] and so we have a formula for the satellite mass: \[M_s = \frac{T(R)R^2}{ \omega_s^2 R^2 - G M_e}\]
Alternatively if we know the mass of the satellite (\(M_s\)) and the distance of the satellite from the centre of the earth \(R\) then we can use equation 5 to evaluate \(T(R)\) and then equation 4 to find the resultant tension at the ground anchor: \(T(r_e)\)

And if we have \(T(r_e)\) and \(M_s\) and we want to evaluate R, then we can use equations 4 and 5 and with the aid of a numerical solver we can evaluate R.

Which ever way we do it, when the tension has been determined, we can return to equation 3 to find the mass per unit length and then integrate to find the total mass of cable required: \[M_c = \int^R_{r_e}dm(r) = \frac{\alpha}{Y} \int^R_{r_e}T(r)dr\] That integral can be done numerically.

I've written a calculator using Google sheets and I've published it as a template. To have a look, log-in using a google (gmail) account and click on this link.

Things to consider:
stability to perturbations including wind
could attached kites be helpful in reducing the maximum tension?
collisions: birds, planes, other satellites, space debris.
safety if it snaps. One option would be for the anchor at the earth to release as soon as a break has been detected.



Notes:
In the calculations above we've used the following:
Gravitational Constant: \(G \approx 6.67 \times 10^{-11} \frac{m^3}{s^2 kg}\)
Mass of earth: \(M_e=5.97 \times 10^{24} kg\)
Radius of earth: \(r_e = 6.37 \times 10^6 m\)

Friday, December 5, 2014

Greek Probabilities

Each month IBM posts up a mathematical puzzle. The November 2014 puzzle was as follows:
There are 24 letters in the Greek alphabet, from alpha, beta, and gamma, through chi, psi, and omega. The names of the letters are of different lengths, from short two-letter names like mu or pi through the seven letters in omicron. We want to generate a random Greek letter, where the probability of a letter is determined by the length of its name, and the following:
Psi is a special case with a probability of zero (we are actually generating one of 23 letters) Letters whose names are composed of an even number of letters have a probability of 2^length times the probability of the letter Tau
All the other letters (except Psi and Tau) have a probability of 2^(length-2) times the probability of Tau So, for example, the letter mu (with a length of two letters) should appear 2^2 = 4 times more than tau. Omicron (with a length of seven letters) should appear 2^(7-2), or 32 times more often than tau.
The task is to generate the random Greek letter using six coin tosses, in which the coins are chosen from two types: one with a probability of p_1 of getting a tails and 1-p_1 of getting heads; and the other with a probability of p_2 for tails and 1-p_2 for heads.


One solution has p_1 = 8 / 17
and p_2 = 1 / 3
We can represent the solution as a circle. We divide up the circumfrence of the circle such that the probability is proportional to the length of the section of the cicumfrence. Then we can present a solution as an image:

Tuesday, August 26, 2014

Pretty pictures with R

One of the nice things about using the R programming language, apart from the fact that it is free, is that it can produce beautiful graphics. As well as do all sorts of standard and non-standard plots and graphs, R can also generate animated GIF files.
For example, suppose we have some height: \[h= exp( - \alpha r ) cos ( \beta \ r + p)\] where r is the distance from the origin: \[r=\sqrt{ x^2 + y^2 }\] Also we have constants
\(\alpha\) which determines the decay
\(\beta\) determines the frequency
and \(p\) is the phase

We generate a matrix which we've called h. Each element corresponds to an x and a y. We can convert that h matrix to an image and then cycle through values of phase to animate it. We will then have an image such as the one above.
Here is the R code that generates the h matrix:
 # generate an image of a ripple  
 rippleMat <- function( numPx, decay, phase, freq){  
   x <- seq( -1, 1, 2/ (numPx-1))  
   y <- x   
     
   # We have a matrix r[i,j] = sqrt( x[i] * x[i] + y[j] * y[j])  
   r <- sqrt( outer( x*x, y*y, "+"))    
     
   h <- exp( - decay * r ) * cos( phase + freq * r )  
   return (h)  
 }  
   
(The full code used to generate that animated gif is show below)

Alternatively, we could also plot a spiral. The basic equation is: \[r = exp( \phi \alpha + p ) \] Where r is the distance from the origin, \(\phi\) is the angle to the x-axis, p is the phase and \(\alpha\) determines how tightly coilded it is.
We could of course plot such a path using a line chart in R. But it might be a bit more interesting to color it in. Again we could iterate through values of the phase p to generate an animated GIF file:

We could also build an RGB color triangle. Each corner is set to one of { red, green, blue }. That color then fades across the triangle until it is completely absent on the opposite edge.

Just in case you're interested, here's the full source code for the triangle:
# Author: Philip Kinlen, Aug 2014
library(grid) # required for grid.raster(.)

# numPix should be an integer,
# and gapX < (numPix / 2) 
# where gapX is the minimum number of horizontal pixels from the triangle to the edge of the image
RGBTriangle <- function(numPix = 100, gapX = 10, doSmoothing=FALSE) {
  # the verticle gap between the triangle and edge of image
  gapY            <- round(numPix * ( 0.5 - sqrt(3)/4) + gapX * sqrt(3)/2)

  xArr            <- 1:numPix
  yArr            <- numPix:1   # the function call grid.raster(..) below will put the element with 
                                # highest row number at the bottom. 
                                # For our array, the last shall be at the bottom.
  # The x's a constant in the matrix columns, y's are constant along the matrix rows.
  xMat            <- matrix(rep(xArr, numPix), numPix, numPix, F)  
  yMat            <- matrix(rep(yArr, numPix), numPix, numPix, T)  

  m1              <- sqrt(3)                        # slope
  c1              <- gapY - m1 * gapX               # intercept
  
  m2              <- -sqrt(3)                       # slope
  c2              <- gapY - m2 * (numPix - gapX)    # intercept
  
  height          <- numPix - 2 * gapY              # height of triangle in pixels

  red             <- matrix(mapply(getRed,   xMat, yMat, m1, c1, m2, c2, gapY, height), numPix, numPix, T)
  green           <- matrix(mapply(getGreen, xMat, yMat, m1, c1, m2, c2, gapY, height), numPix, numPix, T)
  blue            <- matrix(mapply(getBlue,  xMat, yMat, m1, c1, m2, c2, gapY, height), numPix, numPix, T)
  
  col             <- rgb(red, green, blue)
  dim(col)        <- dim(red)
  grid.raster(col, interpolate=doSmoothing) 
  # One possible source of minor confusion is that when we write a cartesian co-ordinate (x,y)
  # the x comes first.
  # but in a matrix, we specify the row first, which corresponds to y
} 

#################################################################
getRed <- function( x, y, m1, c1, m2, c2, gapY, height){
   if ( isInsideTriangle( x, y, m1, c1, m2, c2, gapY)){
     res <- (y - gapY) / height    # A rounding error may cause res to be outside the range [0,1]
     return (max(0, min(1, res)))  # If left untreated that could cause an error in grid.raster(.)
   } else
     return (0)     
}
#################################################################
getBlue <- function( x, y, m1, c1, m2, c2, gapY, height){
  if ( isInsideTriangle( x, y, m1, c1, m2, c2, gapY)){ 
    res <- sqrt(3) * ((y - c2) / m2 - x) / ( 2 * height) 
    return (max(0, min(1, res))) 
  } else
    return (0)     
}  
#################################################################
getGreen <- function( x, y, m1, c1, m2, c2, gapY, height){
  if ( isInsideTriangle( x, y, m1, c1, m2, c2, gapY)){
    res <- 0.5 * (m1 * x + c1 - y) / height
    return (max(0, min(1, res))) 
  } else
    return (0)     
}  
#################################################################
isInsideTriangle <- function(x,y, m1, c1, m2 ,c2, gapY){
  
  if (  ( y >  gapY) 
      & ( y < ( m1 * x + c1 ) ) 
      & ( y < ( m2 * x + c2 ) ) )
    return(T)
  else 
    return (F)
}

And here is the code used to generate the animated GIFs contains the ripples and the spiral
# Author: Philip Kinlen, Aug 2014
###################################################################
library( caTools)  # library required for write.gif(..)

###################################################################
animateRipples   <- function(){
  
  filePath       <- paste(getwd(), "/animatedRipples.gif", sep="")
  numImages      <- 25
  numSecForCycle <- 1
  
  numPx          <- getNumPixels()
  phases         <- -2 * pi * ( 0:(numImages-1)) / numImages
  
  # The following line does the main work generating the data. 
  # The call to sapply(..) does the work that would have been done by a for-loop in
  # a language such as a Java.
  zeds           <- sapply(phases, rippleMat, simplify="array")    

  delayCS        <- 100 * numSecForCycle / numImages  # in hundreths of a second
  
  write.gif( zeds, filePath, col= getCustomColors(), delay = delayCS )
  
}
###################################################################
animateSpiral     <- function(){
  # for MS Windows might need to set the separator argument to a back-slash
  filePath        <- paste(getwd(), "animatedSpiral.gif", sep="/")
  
  numImagesPerSec <- 20
  numSecForCycle  <- 5
  
  numImages       <- numImagesPerSec * numSecForCycle
  numPx           <- getNumPixels()
  
  phases          <- 2 * pi * ( 0:(numImages-1)) / (numImages * getSpiralBeta())

  # The following line does the main work generating the data. 
  # The call to sapply(..) does the work that would have been done by a for-loop in
  # a language such as a Java.
  zeds            <- sapply( phases, spiralMat, simplify="array")
  
  delayCS         <- 100 * numSecForCycle / numImages  # in hundreths of a second
  
  write.gif( zeds, filePath, col= getCustomColors(), delay = delayCS )
}
###################################################################
getSpiralBeta <- function(){
  return ( 0.2)
}
###################################################################
getNumPixels <- function(){
  return ( 200)
}
###################################################################
getCustomColors   <- function(){
  maxCol     <- 255
  
  rising     <- 0:maxCol
  falling    <- maxCol - rising
  
  myRed      <- c(falling, rising) * 0.4 + 0.2
  myGreen    <- c(falling, rising) * 0.4 + 0.2
  myBlue     <- c(falling, rising) * 0.8 + 0.2
  
  return (rgb(  myRed, myGreen, myBlue,  maxCol, maxColorValue=maxCol ))
}  
###################################################################
# generate a matrix representation of a ripple
rippleMat <- function( phase, numPx = getNumPixels(), decay = 2.5, freq = 12 * pi){
   x <- seq( -1,  1, 2/ (numPx-1))
   y <- x 
   
   # We have a matrix r[i,j] = sqrt( x[i] * x[i] + y[j] * y[j])
   r <- sqrt( outer( x*x, y*y, "+"))   
   
   h <- exp( - decay * r ) * cos( phase + freq * r )
   return (h)
}
###################################################################
spiralMat    <- function(phase, beta = getSpiralBeta() , numPx= getNumPixels(), generateImage= FALSE){

  alpha    <- 0.1   # the main spiral follows: r = exp( alpha * phi)
  
  x        <- seq(-1, 1, 2 / (numPx - 1))
  y        <- x
  matR     <- sqrt(outer(x*x, y*y, "+"))
  
  matX     <- matrix(rep(x, numPx), numPx, numPx, F)  
  matY     <- matrix(rep(y, numPx), numPx, numPx, T)  
  theta    <- ( atan2(matY, matX) + pi) + phase

  # Note that:
  #      r >= exp( (2 * pi * (rev-1) + theta) * alpha)  
  # and  r <= exp( (2 * pi *  rev    + theta) * alpha)
  # where rev is the number of revolutions 
  rev      <- ceiling((log(matR)/ alpha  - theta) /( 2 * pi))

  phi      <- rev * 2 * pi + theta
  z        <- sin(phi * beta)

  if ( generateImage)
    image( z,  col = getCustomColors())
  else 
    return (z)
}
###################################################################

Tuesday, July 15, 2014

Colliding Bullets

Each month, IBM publish a mathematical puzzle. For example in May 2014 they posted The Colliding Bullets Puzzle:

Every second, a gun shoots a bullet in the same direction at a random constant speed between 0 and 1.
The speeds of the bullets are independent uniform random variables. Each bullet keeps the exact same speed and when two bullets collide, they are both annihilated.
After shooting 20 bullets, what is the probability that eventually all the bullets will be annihilated?
Please supply the answer rounded to the 10th decimal digit.


It is reasonably straight-forward to let a Monte-Carlo engine work it out.
For example see this link to Java code in github.

However, to get an answer to an accuracy of 10 decimal digits would take a lot of computing power or a more clever algorithm.

Let \(v_i\) be the speed of bullet i
( I’m going to use a zero based numbering system, not because I really like it, but rather to be consistent with the Java code. It makes it simpler to transition between the two. ).

We're assuming that the bullets are infinitesimally small and the the probability of 3 colliding together at the same point in time and space is (almost surely) zero.
Also, we note that collisions can only take place between odd and even bullets. We can't have two even-numbered bullets colliding together.
To explain why, we consider what would happen if two even bullets were to collide. Then we would need all the intermediate bullets to have been annihilated. And the number of bullets between say bullet 2m and 2(m+n) is 2n-1 which is odd. But bullets only annihilate in pairs. So we can't have an odd number of intermediate bullets annihilating. Thus two even bullets cannot collide with each other. Similarly two odd bullets cannot collide with each other.
(Far) below we deal with the case when we have 4 bullets but first we'll look at an more general case, with generalized initial conditions.
Suppose \(2N\) bullets are fired. We will label the starting condition to be \(\mathcal{J}_{2N} (\textbf{X},\textbf{t})\):
each bullet has position \(X_i\) with
\(X_i \ge X_{i+1} \) and \(X_i \ge 0 \ \ \ \ \ \forall i: 0 \le i \le N-1 \)

start times: \(t_i \) with
\(t_i \le t_{i+1} \) and \(t_i \ge 0 \ \ \ \ \forall i: 0 \le i \le N-1 \)

and speeds as before \(v_i \sim U(0,1) \)

We define \(\mathcal{A}_{2N} \) to be the probability that all bullets released in our generalized version will be annihilated.
In the case when there are 2 bullets ( \(N=1 \) )
they will collide when \(v_0 < v_1\).
Otherwise they won't.
So the probability that they annihilate each other is: \[\mathcal{A}_{2} = \frac{1}{2}\] Now consider the case when we have an even number of bullets and there are more than 2 of them, i.e \(N>1\):
The only case when there are no collisions at all is when the velocities are strictly ordered: \[v_i > v_{i+1} \ \ \ \ \ \forall i: 0 \le i < N-1 \] In that case, the first bullet is the fastest and the last is the slowest.
On the other hand, if last is not the slowest, then we know there will be at least one collision.
The probability that the last is not the slowest is: \(p_{2N} = \frac{2N-1}{2N} \).
In this case, when the last is not the slowest there will a first collision. We don't know which pair will collide first, but we know that the first collision will definitely occur. After this first collision occurs, we then have 2 fewer remaining bullets.
Immediately after the first collision we will then have a valid starting condition \(\mathcal{J}_{2N-2}(\textbf{X}',\textbf{t}') \). If all the original 2N bullets are to be annihilated, then we will need the remaining (2N-2) bullets to be annihilated and so we find: \[\mathcal{A}_{2N} = \frac{2N-1}{2N} \mathcal{A}_{2N-2}\] Using this recursive definition we find: \[\mathcal{A}_{2N}=\prod_{n=1}^N\frac{2n-1}{2n}\]

Alternative approaches

Integration

Another way to try to get the solution is to analytically do integrals over the probability density functions and thus end up with an exact solution. Now lets move on to the case when we have 4 bullets.
We'll introduce some notation. The probability that all 4 bullets will be annihilated is:
\[A_4= \Phi (**,|/) + \Phi (\bigwedge,\bigwedge) + \Phi (|/,\bigwedge) \] Where:
\(\Phi (**,|/)\)is the contribution from the case when we don’t have any restrictions on the first pair and \(v_2 < v_3\).
I.e. the second pair is NOT on a collision course.

\(\Phi (\bigwedge,\bigwedge)\) is the contribution from the case when \(v_0 < v_1\) and \(v_2 < v_3\).
I.e. that's the case when the first are on a collision course and so are the second pair. Though it doesn't necessarily mean that the first pair will collide with each other.

\(\Phi (|/,\bigwedge)\) is the contribution from the case when \(v_1 < v_0\) and \(v_2 < v_3\).
i.e. the first pair are NOT on a collision course, but the second pair are.

First lets evaluate \(\Phi (**,|/)\). In this case the second pair are not on a collision course.
For all to be annihilated we must have a collision between bullets 1 and 2 and then a collision between 0 and 3.

We first work out the minimum velocity of \(v_0\) required to ensure that bullets 0 and 1 do not collide before the collision of bullets 1 and 2.
With a little bit of algebra ( and or geometry ) we find that minimum velocity to be: \[\beta=\frac{v_1 v_2}{2 v_2 - v_1}\] So we have the constraints:
\[0 < v_2 < 1 \] \[0 < v_1 < v_2 \] \[\beta < v_0 < v_2 \] \[v_0 < v_3 < v_2 \]
Note that the probability density function of each of the velocities \(v_i\) is 1 inside the range 0 to 1. And it is zero outside.
The (subtle) difference between the condition \( v_2 < 1 \) and \( v_2 \leq 1 \) is on a set of zero measure which won't affect our results.
Integrating over the probability space we find: \[\Phi (**,|/)=\int_0^1 dv_2 \int_0^{v_2}dv_1 \int_{\beta}^{v_2}dv_0 \int_{v_0}^{v_2}dv_3 \] Doing the \(dv_3\) integral we find: \[\Phi (**,|/)=\int_0^1 dv_2 \int_0^{v_2}dv_1 \int_{\beta}^{v_2}dv_0 ( v_2 - v_0) \] No doing the \(dv_0\) integral we find: \[\Phi (**,|/)=\int_0^1 dv_2 \int_0^{v_2^2}dv_1 \left( \frac{v_2}{2} - \frac{v_1 v_2^2}{2v_2 - v_1} + \frac{1}{2} \left( \frac{v_1 v_2}{2v_2 - v_1} \right)^2\right) \] To evaluate the \(dv_1\) integral takes a bit of effort.
After doing the work we find: \[\Phi (**,|/)=\int_0^1 dv_2 ( v_2^3 [ 3 - 4 ln(2)]) \] and so: \[\Phi (**,|/)= \frac{3}{4} - ln(2) \approx 0.0568528 \]
Now we’ll evaluate \(\Phi (\bigwedge,\bigwedge)\).
In this case \(v_0 < v_1\) and \(v_2 < v_3\).
Suppose it is also true that \(v_1 < v_2\).
In that case we can deduce that \(v_0 < v_3\)
Thus if \(v_1\) and \(v_2\) collide then so too will \(v_0\) and \(v_3\).
And so all 4 will have been annihilated.

On the other hand if \(v_1\) and \(v_2\) don't collide, then we know that the first pair will collide and so will the second pair.
Either way all 4 bullets will be annihilated.
\[Prob(v_0 < v_1)= \frac{1}{2} \] and \[Prob(v_2 < v_3)= \frac{1}{2} \] So \[\Phi (\bigwedge,\bigwedge)=\frac{1}{2} \frac{1}{2} = \frac{1}{4}\]
Now we move on to: \(\Phi (|/,\bigwedge)\)
We're dealing with the case when the first pair is not on a collision course, but the second pair are. For all 4 bullets to be annihilated we requrie bullets 1 and 2 to collide and then bullets 0 and 3 to collide. We have conditions: \[v_1 < v_0 < v_3 \] \[v_1 < v_2 < v_3 \] However, if we just integrate over those ranges, then we'll have over counted. Because we will have also included a contribution from the case when bullets 2 and 3 collide before bullets 1 and 2. And when that happens bullets 0 and 1 won't be annihilated. So we need to subtract a bit which we'll label \(\phi\)
Thus: \[\Phi (|/,\bigwedge) =\left( \int_0^1 dv_3 \int_0^{v_3}dv_1 \int_{v_1}^{v_3}dv_0 \int_{v_1}^{v_3}dv_2 \right) - \phi \] doing that quadruple integral we find: \[\Phi (|/,\bigwedge) =\frac{1}{12} - \phi \]


Showing the over-counting when bullets 0 and 1 are not annihilated.

Now we turn to that extra bit ( \(\phi\) ).
We have the same constraints as above: \[v_1 < v_0 < v_3 \] \[v_1 < v_2 < v_3 \] but we have some extra constraints, which we'll now look into. We work out the time of collision between 1 and 2 to be: \[T_{12}=\frac{2v_2 - v_1}{v_2 - v_1}\] In this extra bit, that we've over-counted, bullet 3 will have already overtaken bullet 2 before the collision between 1 and 2. I.e. bullet 3 will have passed the (1,2) collision at time \(T_{12}\):
\[v_2(T_{12} -2) < v_3 ( T_{12} - 3 )\] That implies: \[ v_1 v_2 < v_3 ( 2v_1 - v_2) \] and we must have: \[2 v_1 > v_2 \] and \[v_3 > \frac{v_1 v_2}{2v_1 - v_2}\]
The lower limit of the integral over \(v_3\) will be \(\frac{v_1 v_2}{2v_1 - v_2} \), so we need to ensure that: \[ \frac{v_1 v_2}{2v_1 - v_2} < 1 \] and so we require: \[v_1 > \frac{v_2}{2-v_2} \]
Thus \[\phi = \int_0^1 dv_2 \int_{\frac{v_2}{2-v_2}}^{v_2}dv_1 \int_{\frac{v_1 v_2}{2v_1 - v_2}}^{1}dv_3 \int_{v_1}^{v_3}dv_0 \] So: \[\phi = \int_0^1 dv_2 \int_{\frac{v_2}{2-v_2}}^{v_2}dv_1 \int_{\frac{v_1 v_2}{2v_1 - v_2}}^{1}dv_3 (v_3 - v_1) \] and hence: \[\phi = \int_0^1 dv_2 \int_{\frac{v_2}{2-v_2}}^{v_2}dv_1 \left( \frac{1}{2} - v_1 + \frac{v^2_1 v_2}{2v_1 - v_2} - \frac{1}{2} \left( \frac{v_1 v_2}{2v_1 - v_2} \right)^2 \right) \] It does indeed take a bit of work to show that it evaluates to: \[\phi = \int_0^1 dv_2 \left( \frac{v_2}{2} - \frac{5v_2^2}{8} + \frac{v_2^3}{2} + \frac{v_2^3-4v_2}{8(2-v_2)} + \frac{2v_2^2-3v_2^3+v_2^4}{4 (2-v_2)^2} \right)\] (To work out the integral this link may be helpful. )
We find: \[\phi= \frac{6}{24} - \frac{5}{24} + \frac{3}{24} - \frac{4}{24} + \frac{17}{24} - ln(2)\] So we end up with: \[\phi= \frac{17}{24} - ln(2)\] and we now have: \[ \Phi (|/,\bigwedge) = \frac{1}{12} - \left( \frac{17}{24} - ln(2) \right) = ln(2) - \frac{5}{8} \] We return to the probability that all 4 bullets will be annihilated: \[A_4= \left[ \Phi (**,|/) \right] + \Phi (\bigwedge,\bigwedge) + \left\{ \Phi (|/,\bigwedge) \right\} \] We now have evaluated all the terms on the right hand side: \[A_4= \left[ \frac{3}{4} - ln(2) \right] + \frac{1}{4} + \left\{ ln(2) - \frac{5}{8} \right\} \] And so: \[A_4 = \frac{3}{8}\]
So we now know that, when 4 bullets are fired, the probability that all will be annihilated is \(\frac{3}{8}\)

Let's not integrate:

I am aware that for some people, when they are imagining a fun evening, thoughts of doing quadruple integrals by the fireside do not enter their mind. For such individuals, we offer an alternative approach to working out \(A_4\) without doing any integration at all.
We write the probability that all 4 will be annihilated: \[A_4=\Phi \left( \bigwedge , \bigwedge \right) + \Phi \left( |/ , |/ \right) + \Phi \left( |/ , \bigwedge \right) + \Phi \left( \bigwedge, |/ \right) \]
\( \Phi \left( \bigwedge , \bigwedge \right) \) is the contribution from the case when
\( (v_0 < v_1 ) \) and \( (v_2 < v_3 ) \)
We have already worked that out above, without using integration and found it to be: \[ \Phi \left( \bigwedge , \bigwedge \right) = \frac{1}{4} \]
Now moving on to \( \Phi \left( |/ , |/ \right)\) which is the contribution from the case when
\( (v_1 < v_0 ) \) and \( (v_3 < v_2 ) \)
When that is true, all 4 bullets will be annihilated when we have the additional conditions:
\( (v_0 < v_3 ) \) and \( (v_1 < v_3 ) \)
So we have the conditions: \[ v_1 < v_0 < v_3 < v_2 \] There are \(4!=24\) ways in which the 4 \(v_i \)s can be ordered. Each has the same probability \( \frac{1}{24} \). Only 1 of the 24 permuations obeys the conditions that we have here and so we can deduce that: \[\Phi \left( |/ , |/ \right) = \frac{1}{24} \]
Now moving on to: \( \Phi \left( |/ , \bigwedge \right) \). Remember, it is the contribution to \(A_4 \) from the case when the first pair not on a collision course, but the second pair are. So we have conditions:
\( (v_1 < v_0 ) \) and \( (v_2 < v_3 ) \)
for all to be annihilated, we have additional constraints: \( (v_0 < v_3 ) \) and \( (v_1 < v_2 ) \)
and we also require that bullet 2 does not collide with 3 before 1.
After a bit of work, we find that that condition is: \[ 2 v_1 v_3 < v_2 ( v_1 + v_3) \] We'll denote the probability to be \(P\left\{.\right\}\)
So we have: \[ \Phi \left( |/ , \bigwedge \right) = P \left\{ (v_1 < v_0 ) \& (v_2 < v_3 ) \& (v_0 < v_3 ) \& (v_1 < v_2 ) \& \left[ 2 v_1 v_3 < v_2 ( v_1 + v_3) \right] \right\} \] We'll return to that in a moment. But first we'll examine \( \Phi \left( \bigwedge, |/ \right) \). Using a similar argument we find that: \[ \Phi \left( \bigwedge, |/ \right) = P \left\{ (v_0 < v_1 ) \& (v_3 < v_2 ) \& (v_0 < v_3 ) \& (v_1 < v_2 ) \& \left[ 2 v_0 v_2 > v_1 ( v_0 + v_2) \right] \right\} \] We now note that the \(v_i\) are independent, identically distributed random variables and so we could permute the indices on the right hand side of the above equations without changing the value. On the right hand side of the equation immediately above change: \[i_{new} = (i_{old} + 1) mod \ 4\] And so we have: \[ \Phi \left( \bigwedge, |/ \right) = P \left\{ (v_1 < v_2 ) \& (v_0 < v_3 ) \& (v_1 < v_0 ) \& (v_2 < v_3 ) \& \left[ 2 v_1 v_3 > v_2 ( v_1 + v_3) \right] \right\} \] Now note that the right hand side that equation contains a term: \[ 2 v_1 v_3 > v_2 ( v_1 + v_3)\] And the formula above for \( \Phi \left( |/, \bigwedge \right) \) contains the same term with the '>' switched to a '<'. When we're working out the probabilities we can ignore '=' case, since that is on a set of zero measure.
Now note that generally speaking: \[P\left\{K \& L \right\} + P\left\{K \& \ not(L) \right\} = P\left\{K \right\}\] And so we find: \[\Phi \left( \bigwedge, |/ \right) + \Phi \left( |/, \bigwedge \right) = P \left\{ (v_1 < v_0 ) \& (v_2 < v_3 ) \& (v_0 < v_3 ) \& (v_1 < v_2 ) \right\}\] Of the 24 permutations of \(v_i\), we find that just 2 obey those conditions and so \[\Phi \left( \bigwedge, |/ \right) + \Phi \left( |/, \bigwedge \right) = \frac{2}{24} = \frac{1}{12}\] Note that we have not worked out what each of the terms on the left hand side is equal to, but we have worked out their sum, which is what we needed.
Now returning to \(A_4\) we find: \[A_4=\frac{1}{4}+ \frac{1}{24}+ \frac{1}{12} = \frac{3}{8}\]