<< Algorithms
Random Yet Smooth Motion

A point moves about the plane. We seek a motion of the point that is both random and at the same time smooth. Those characteristics would seem to be contradictory. If the motion is smooth then given a position of the point, its position a little later would be in nearly the same direction as the point had been traveling. A point’s immediate future positions would be much more predictable than if the motion were not smooth.

Here is how to do it. First consider a random walk on the line, where a step is taken in one direction or the other at random. Place the line vertically so its trajectory can be graphed with time as the horizontal axis. In the following diagram the graph of the random walk,  X(t),  is colored red. Then take the running average over a number of steps, say 10. Better, take the weighted average using a bell-shaped function. In the diagram the graph of the resulting smooth motion  X̅(t)  is colored blue.  In signal jargon, the running average is a kind of low pass filter. Mathematically the running average is the convolution of the random-walk with the weight function.

For the plane, do such a convolution on each dimension separately. Start with the “drunkard’s walk”:  choose a random direction and move one step in that direction, then repeat. Now consider the projection of that motion onto the X and Y axes, so you have two functions of time  X(t), Y(t).  Consider each function separately. The graph is angular and bumpy. After the walk has progressed for a certain number of steps – say 10 – start taking the running average weighted by a bell shaped curve whose width covers 10 steps (in other words, convolve each of  X(t)  and  Y(t)  with a bell shaped curve 10 steps wide). That smoothes the two functions. Then combine them again to trace out a random but smooth motion   X̅(t), Y̅(t).

The moving point might slow, come to a stop, reverse direction and accelerate. In that case, though the path is not smooth at the turning point the motion is. The qualification smooth applied to the motion of a point means that its velocity varies continuously.

A suitable weight function is positive within the interval 0 to 1, zero outside it, and its integral is one. In use the interval can be scaled to cover the number of steps required. Example weight functions:

The first example is the cosine function in the interval  –π to π,  suitably scaled and shifted:
      1 + cos  (x  ½)
for  x  between  0  and  1.

The second example is
      A e 1 / ( (2x – 1) 2 – 1)
where A is about 4.5 so the integral from 0 to 1 is 1. When  x = ½  this function is about 1.66 so its graph is a bit fatter than that of the first weight function above. By replacing the exponent 2 with a larger, even, exponent, you get a still fatter graph.

The third is the “rectangular step function.” It gives the uniformly weighted average.

It turns out that any more or less bell shaped curve will do for the weight function. Even the rectangular function, in other words an ordinary average instead of a weighted one, makes little difference in the general look of the resulting curve.  The program below uses the first weight function graphed above.

A novel feature of the algorithm, found among the papers of M. Ralston, presented in the program below is that the smooth trajectory is computed as the random walk progresses.  After computing an initial  ‘F’  steps, the smooth trajectory is prolonged one step at a time.

Because we choose a random direction for each step of the random walk and then project onto X and Y, the projected length of the step varies, between zero and the step length in the plane. The average projected length is  2 / π  times the planar step length.

In the program below, the moving point is prevented from wandering outside the graphic window as follows. If a random step would go outside the window, another random step is tried until it doesn’t. (Without this restriction, in dimensions one or two a random walk eventually returns arbitrarily close to its starting point, and every other point in the plane. In three or higher dimensions it might wander off never to return.)

The source code can be modified for Visual Basic, FORTRAN, or any other BASIC-like compiler.  An executable (for Windows computers) with an added interface that allows you to change various parameters without re-compiling, can be downloaded from here.

If we choose for the weight function the rectangular function, that is, we take the simple running average, then  bumpf()  below is identically one. In that case its declaration, preliminary computation, and all  “* bumpf(...)”  can be removed.

' Defining constants ...
' F is the most important. STEPLEN should be smaller the larger F is, to
' better fit in the window.  M should be about STEPLEN for a solid line.
Const F = 3                           'number of steps to average over
Const STEPLEN = 150                   'step length
Const M = 150                         'number of points in one step of the random walk
' larger F, fewer tight loops;
' larger STEPLEN, magnify drawing (also points of trajectory further apart);
' larger M, points of trajectory closer together;

' Derived constant
Const N = M * F                       'total number of points to average over

' Speed constant
Const DELAY = 2                       'the larger this is the slower the random walk

' Compiler switch
#Const showsteps = -1                 'show the random walk or not

Function PBMain
 Local i, j, k As Long                'indices for the most part
 Local winC As Dword                  'handle of graphic window
 Local xbrown0, ybrown0 As Single     'beginning of brownian step
 Local xbrown1, ybrown1 As Single     'end
 Local limit As Long                  'half screen height
 Local x, y As Single                 'current point
 Local dx, dy As Single               'change in direction of step
 Local th As Single                   'angle of a step
 Local pd As Single                   'pi doubled
 Dim xline(M - 1) As Single, yline(M - 1) As Single    
 Dim xbrown(N - 1) As Single, ybrown(N - 1) As Single  
 Dim bumpf(N - 1) As Single           'weight function

 ' Set up graphic window, a square window that fills the screen vertically
 Desktop Get Client To i, j             'get computer’s screen size, horizontal & vertical
 Console Set Loc 0, j + 1               'hide console window
 Graphic Window "", (i - j) / 2,0, j,j To winC    'create graphic window
 Graphic Attach winC, 0, ReDraw         'choose it, specify drawing not shown until ReDraw
 limit = j / 2                          'half the window height
 Graphic Scale (-limit, limit) - (limit, -limit)  '(0, 0) at center of graphic window

 ' Instructions to user
 Graphic Set Focus
 Sleep 10
 Graphic Font "calibri", 15, 0
 Graphic Print
 Graphic Print "       To start a new path, press Spacebar."
 Graphic Print
 Graphic Print "       Now press Spacebar, or Esc to quit."
 Graphic ReDraw
 Do                                     'place and keep console in focus
  Console Set Focus
  Sleep 100
 Loop Until InStat
 If WaitKey$ = $Esc Then Exit Function  'wait for a key to be pressed

 Randomize Timer                        'random seed for Rnd

 pd =  8 * Atn(1)                       'compute 2 * pi

 ' Pre-compute the weight function.  Note that  (i - N / 2) / N   =   i / N - 1/2.
 ' Because the integral from 0 to 1 of the function  1 + Cos 2(pi)(x - 1/2)
 ' equals 1, the sum of  bumpf(0 to N - 1)  equals, to many decimal places, N.
 ' In use, after multiplying by the weight function and summing, divide by N.
 For i = 0 To N - 1
  bumpf(i) =  1 + Cos(pd * (i - N / 2) / N)

 Do                                     'new path

  xbrown0 = 0 :ybrown0 = 0
  Reset xbrown(), ybrown()
  Graphic Clear

  Graphic Ellipse (-2, 2) - (3, -3), %Red, %Red           'draw starting dot
  Graphic ReDraw
  Sleep 200                                               'slight pause

  Do                                                      'F random walk steps

    th = pd * Rnd                                         'between 0 and two pi
    dx = Cos(th) * STEPLEN                                'random direction,
    dy = Sin(th) * STEPLEN                                'keep trying until within window
    xbrown1 = xbrown0 + dx
    ybrown1 = ybrown0 + dy
   Loop Until xbrown1 > -limit And xbrown1 < limit And ybrown1 > -limit And ybrown1 < limit

   For j = 0 To M – 1                                     'step from brown0 to brown1
    xline(j) = xbrown0 + dx * j / M
    yline(j) = ybrown0 + dy * j / M

   #If showsteps Then                                     'show random walk in light gray
   Graphic Line (xbrown0, ybrown0) - (xbrown1, ybrown1), &Ha0a0a0
   Graphic ReDraw
   #End If

   If DELAY Then Sleep DELAY

   xbrown0 = xbrown1
   ybrown0 = ybrown1

   ' Compute running average over section, incremented over one piece.
   For k = 0 To M - 1
    x = 0
    y = 0
    For i = 0 To k - 1
     x = x + xline(i) * bumpf(i + N - k)
     y = y + yline(i) * bumpf(i + N - k)
    For i = k To N - 1
     x = x + xbrown(i) * bumpf(i - k)
     y = y + ybrown(i) * bumpf(i - k)

    'draw point
    Graphic Set Pixel (x / N, y / N), %Blue
    Graphic ReDraw

   ' Discard first step of brown() and tack on line() step.
   For i = M To N - 1                  'move back by M
    xbrown(i - M) = xbrown(i)          'brown(0 to N - M - 1) = brown(M to N - 1)
    ybrown(i - M) = ybrown(i)
   For i = N To M + N - 1
    xbrown(i - M) = xline(i - N)       'brown(N - M to N - 1) = line(0 to M - 1)
    ybrown(i - M) = yline(i - N)

   Console Set Focus                   'in case user clicked on the graphic window

  Loop Until InStat                    'stop if key pressed

 Loop Until WaitKey$ = $Esc            'new path
End Function

A path in three dimensions can be generated similarly. For the random step direction, instead of choosing a direction indicated by a random angle, which amounts to a point on the unit circle, choose a random point on the unit sphere.

The direction must be uniformly distributed; thus a pair of latitude and longitude each chosen randomly separately cannot be used (the result would concentrate toward the poles). An amazing geometrical fact discovered by Archimedes, the basis of what is now called the Hatbox Theorem, makes the solution easy.

Imagine a cylinder circumscribing the unit sphere, aligned with the vertical z axis. Then given any patch of area on the sphere, its projection onto the cylinder – by horizontal lines from the z axis – has the same area. (In cartography this is known as Lambert’s cylindrical projection.) Thus a uniform powder of points on the sphere projects to a uniform powder of points on the cylinder. Of course that is true in reverse, a uniform powder on the cylinder projects to a uniform powder on the sphere.

The cylinder can be slit vertically and unrolled without distortion into a 2 by 2π rectangle. A point on this rectangle can be chosen randomly and uniformly by so choosing its rectilinear coordinates separately.

Here’s how to modify the above program for three dimensions. Declare new variables  z  after  y,  and after  ybrown0,  ybrown1,  ybrown(),  and  yline()  then duplicate the appropriate instructions. To obtain a random point uniformly distributed over the unit sphere, replace the code which uses one call to Rnd:
    th = pd * Rnd
    dx = Cos(th) * STEPLEN                                'random direction,
    dy = Sin(th) * STEPLEN                                'keep trying until within window
    xbrown1 = xbrown0 + dx
    ybrown1 = ybrown0 + dy
   Loop Until xbrown1 > -limit And xbrown1 < limit And ybrown1 > -limit And ybrown1 < limit
with this code which uses two:
 Local dx, dy, dz As Single      'direction of step
 Local s As Single               'temporary variable
    dz = Rnd * 2 - 1             'between -1 and 1
    s = Sqr(1 - dz * dz)         'radius of slice at dz
    th = pd * Rnd                'between 0 and two pi
    dx =  s * Cos(th)
    dy =  s * Sin(th)
    xbrown1 = xbrown0 + dx * STEPLEN
    ybrown1 = ybrown0 + dy * STEPLEN
    zbrown1 = zbrown0 + dz * STEPLEN
   Loop Until xbrown1 > -limit And xbrown1 < limit And ybrown1 > -limit And ybrown1 < limit  And zbrown1 > -limit And zbrown1 < limit