The Magnetic Pendulum

Demonstrating the butterfly effect with a magnetic pendulum

Introducing the Magnetic Pendulum

This article is based on an article named "Experimente zum Chaos" (translated.: "Experiments on Chaos") [1] from the the German magazine Specktrum der Wissenschaft (German version of the Scientific American magazine). The article dates back to 1994, and introduces among other things a simulation showing the chaotic motion of a pendulum under the influence of gravity and three magnets. For some reason, this simulation fascinated me and so I wrote a program implementing it. This article will describe the program and the science behind simulating the magnetic pendulum.

Video 1: The animation shows color patterns created by a simulation with 6 rotating magnets.

The simulation is based on computing the motion of a metallic pendulum under the influence of three magnets. A color image is obtained by integrating the equations of motion for all possible starting points in a two dimensional grid and coloring the starting point with the color of the magnet over which the pendulum came to a rest. Since this involves lengthy integration of the pendulum trail it may take some time to complete the calculation. A few hours is not unusual for a typical image in HD size (1980 x 1080 pixels). On a modern GPU such a task could be done in real time but this simulation does not make use of this technology.

What you need to enjoy and understand this article: What you get:

The Butterfly Effect

The simulation is a good example for the so called butterfly effect. For those of you that are unfamiliar with that phrase, here is the brief explanation taken directly from Wikipedia [1]:

The butterfly effect is a phrase that encapsulates the more technical notion of sensitive dependence on initial conditions in chaos theory. Small variations of the initial condition of a dynamical system may produce large variations in the long term behavior of the system.

So how does this all apply to the magnets and pendulum simulation? We will come to that later. We will examine the equations of motion, the integration algorithms and the software itself. But before we do I will show you something even better: A video explaining it indepth. The following video was created by Aldo Cavini Benedetti who did a great job in explaining what the magnets and pendulum simulation is all about. If you watch to the end you will will see that it's actually so simple that even cats understand it:

Video 2: Chaos Theory explained; © 2007 Aldo Cavini Benedetti

Simulation overview

By now you should have an idea where we are heading. The next sections describe a numerical model that demonstrates how small changes in the initial conditions of the simulation can result in large variations of the results. The result is an unpredictability of the simulation result since even the smallest change in the environment might effect the outcome dramatically.

Let's start with the details. The classical model assumes having a magnetic pendulum which is attracted by three magnets with each magnet having a distinct color. The magnets are located underneath the pendulum on a circle centered at the pendulum mount-point. They are strong enough to attract the pendulum in a way that it will not come to rest in the center position. The following picture illustrates the model when watched from above. Colored circles symbolize magnets, the cross in the middle the pendulum mount point. The simulation will calculate the route taken by the pendulum under the influence of all three sources plus gravity and friction. Due to energy loss caused by friction, the pendulum will earlier or later stop over one of the magnets. The starting point is then colored with the color of this very same magnet. Doing this for all pixels will result in a pretty map showing a pattern composed of red, green, and blue pixels.

Pendulum traces above the magnets
Image 1 (left): Rendering of the experimental setup (courtesy of Paul Nylander [3]).
Image 1 (right): Small variations in the initial positions lead to large variations in the pendulum movement.

The results are shown in the images 1 and 3. The white spot in the middle of image 3 is caused by the pendulum coming to rest in the middle. This is due to the parameters used in this example. Since I'm not a big fan of low color images, adding more color seemed to be necessary. The only information available for each pixel is the color of the magnet, right? And, if you search the internet, you will quickly find that most codes implementing this simulation stop here, but a little more code can greatly improve the results. Remember, for each pixel, we are calculating the whole trace of the pendulum. So it may be a good idea to translate some of the information from the trace into color information. The most obvious is of course the length of the trace. So using the trace length for determining the pixel brightness seems a logical step. (OK, I admit this is not my idea since it was already done in the original version published in Scientific American.) Doing this can be done by using color scaling functions, taking both trace length and maximum trace length as the parameters. The scaling is applied to the source color. In general, the application can use any function, but according to my experience, the following three are useful:

Color functions
Image 2: Overview over different color functions.

As you can see, the result is quite interesting. At least if you have an interest in chaos theory. (If not I wonder why you are still reading.) Starting points resulting in longer traces are shown in darker colors, adding additional complexity to the image. The application allows you to define custom formulas for the color scaling since they will be interpreted using a math expression parser.

Simulation setup overview
Image 3 (left to right): Simulation setup; Color determined by magnet index; Colored determined by magnet index and trace length


The pendulum is a simplified 2D version, assuming the force pulling it back to the centre is following Hookes law (proportional to the distance). This is a simplification, sparing me the effort of calculating rotation angles, cross products, and the whole stuff I would need otherwise. If you won't tell anyone, I could tell you that implementing it physically correct would not be that much additional work, and in fact, I once made a version doing this. But don't forget, my primary objective was getting a picture for my wall, and the physically correct version would have to be mapped to a sphere not a plane. Since I can't hang a sphere on my wall, I'll stick to the 2D version. Of course, the 2D version is valid for small elongations only.

Force distance relationships
Image 4: Example curves of Force vs. Distance for Hooke's and Magnetic Forces.

Magnets are assumed to cause a force proportional to the inverse square of the pendulum distance. In principle, this is akin to the Law of gravity or Coloumbs law. All those laws are very similar, but of course, here we are dealing with (hypothetic) magnetic monopoles, not masses or charges. That assumption is in line with what everyone does when it comes to the pendulum and magnets simulation. In reality, Magnets are dipoles. A dipole causes forces proportional to 1/r³ rather than 1/r². The force calculation does not take this into account although simulating a dipole by two monopole sources would be an option too. The Pendulum is assumed to be made up of iron neglecting eddy currents that would be induced in reality.

Equations of motion

The pendulum movement is calculated by integrating twice over the accelerations acting on the pendulum. Normally, one would not talk about accelerations but forces. According to Newton's second law of motion, the force necessary to move a body equals mass times acceleration. We solve that equation for the acceleration:

\begin{equation} \nonumber\vec{F} = m \vec{a} \end{equation} \begin{equation} \nonumber\vec{a} = {\vec{F} \over m} \end{equation}

Since our initial conditions provide a starting position and a starting velocity (assumed to be null), all we need is to calculate accelerations. For simplicity, mass is assumed to equal one mass unit. Talking about units, I should mention that the simulation in general does not care much about physical units. This is no problem since using real units would just impose scaling factors on the parameters. The following equations list all accelerations relevant for the simulation:

Acceleration caused by the pullback of the pendulum:

\begin{equation} \nonumber\vec{a}_g = k_g \cdot \vec{r} \end{equation}

Accelleration caused by a single magnet:

\begin{equation} \nonumber\vec{a}_m = k_m \cdot \frac{\vec{r}}{|\vec{r}|^3} \end{equation}

Decelleration caused by friction:

\begin{equation} \nonumber \vec{a}_f = -k_f \cdot \vec{v} \end{equation}

Total acceleration of the pendulum:

\begin{equation} \nonumber \vec{a}_t = \vec{a}_g + \Big( \sum_{m=1}^3{\vec{a}_m} \Big) - \vec{a}_f \end{equation}


Integration of the Pendulum trace

Given an initial pendulum position and an initial pendulum velocity, all that is left to do is find a suitable integration scheme and follow the pendulum's trail. For this simulation, the Beeman integration algorithm was used. Applying this scheme does not require much code, and it is pretty accurate. Putting the algorithm into pseudo code looks like this:

while tracing pendulum
    position += Velocity
    acceleration = 0

    for all force sources
    acceleration += acceleration_caused_by_source

    if (pendulum is close to source and velocity is small) then
      stop_magnet = source_index
    end if
    end for

    acceleration -= acceleration_caused_by_friction
    velocity     += acceleration
    trace_length += length(velocity)

    store stop_magnet
    store trace_length
  end while

Implementing the code for tracing the route taken by the pendulum into C++ looks like:

for (int ct=0; ct<m_nMaxSteps && bRunning; ++ct)
  // compute new position
  pos[0] += vel[0]*dt + sqr(dt) * (2.0/3.0 * 
           (*acc)[0] - 1.0 / 6.0 * (*acc_p)[0]);
  pos[1] += vel[1]*dt + sqr(dt) * (2.0/3.0 * 
           (*acc)[1] - 1.0 / 6.0 * (*acc_p)[1]);

  (*acc_n) = 0.0;    // reset accelleration

  // Calculate Force, we deal with Forces proportional
  // to the distance or the inverse square of the distance
  for (std::size_t i=0; i<src_num; ++i)
    const Source &src( m_vSources[i] );
    r = pos - src.pos;
    if (src.type==Source::EType::tpLIN)
      // Hooke's law:          _
      //         _             r         _   
      //  m * a = - k * |r| * --- = -k * r
      //                      |r|
      (*acc_n)[0] -= src.mult * r[0];
      (*acc_n)[1] -= src.mult * r[1];
      // Magnet Forces: _
      //      _         r
      //  m * a = k * -----
      //               |r³|
      double dist( sqrt( sqr(src.pos[0] - pos[0]) +  
                         sqr(src.pos[1] - pos[1]) + 
                         sqr(m_fHeight) ) );
      (*acc_n)[0] -= (src.mult / (dist*dist*dist)) * r[0]; 
      (*acc_n)[1] -= (src.mult / (dist*dist*dist)) * r[1]; 

    // Check abort condition
    if (ct>m_nMinSteps && abs(r)<src.size && abs(vel)<m_fAbortVel)
      bRunning = false;
      stop_mag = (int)i;
  }   // for all sources

  // 3.) Friction proportional to velocity
  (*acc_n)[0] -= vel[0] * m_fFriction;
  (*acc_n)[1] -= vel[1] * m_fFriction;    

  // 4.) Beeman integration for velocities
  vel[0] += dt * ( 1.0/3.0 * (*acc_n)[0] + 5.0/6.0 * 
                 (*acc)[0] - 1.0/6.0 * (*acc_p)[0] );
  vel[1] += dt * ( 1.0/3.0 * (*acc_n)[1] + 5.0/6.0 * 
                 (*acc)[1] - 1.0/6.0 * (*acc_p)[1] );

  // 5.) flip the acc buffer
  tmp = acc_p;
  acc_p = acc;
  acc   = acc_n;
  acc_n = tmp;

You might also like: