Home
Home
Deutsch
German
beltoforion.de
 Contact |  Privacy policy |  CSS |  XHTML

Contents

Download

The magnetic pendulum fractal

Demonstrating the butterfly effect with a magnetic pendulum

Introduction

A few years ago, my physics prof gave me an article named "Experimente zum Chaos" (translated.: "Experiments on Chaos") from the the German version of the Scientific American. 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. The original program was coded quite some time ago, and since then I almost forgot about it. But because now I wanted a high resolution image of this fractal for my wall, I had to recode it. So here it is...

I should mention that the calculation takes quite some time, approx. 4-5 hours is not unusual for an image size of 1000 x 1000 pixels using a moderately fast single core processor. So for you people hoping to see an application that provides quick results in real-time, sorry I can't deliver that! Maybe I should give you a brief overview over what you can expect to find here:

What you need: What you get: Since this article war originally published on CodeProject, so you have some reason to believe that you can find code here that you can actually reuse in your own projects. If you are looking for this, I could tell you that I think that the utility classes used in this project can be quite useful. There is a string tokenizer and other useful utility functions, there is a fixed size vector class using template meta-programming techniques, and there is a class for reading INI files. Pretty useful stuff, but I can't provide you with more details here. This is just not part of this article.

The Butterfly Effect

The program demonstrates the butterfly effect. For those of you that are unfamiliar with that phrase, here is the brief explanation taken directly from Wikipedia:

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? Ok, I'll will show you how, i will give you all the equations the explanations, the algorithms and the software. But before i 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: 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 (image by and courtesy of Paul Nylander).
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

Simplifications

If you ask yourself how I'm implementing all these three dimensional physics, I can tell you it is simple, it is so simple that it is almost embarrassing to tell you. I'm cheating!

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 First Law of Motion, the force necessary to move a body equals mass times acceleration. We solve that equation for the acceleration.

Newtons first law

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:

Equations of motion

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
      break
    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];
    } 
    else
    {     
      //---------------------------------------
      // 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;
      break;
    }
  }   // 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;
}

Using the application

This application needs the Visual Studio 2005 redistributables in order to run properly. Explaining the GUI is simple. There is none! All simulation parameters are read from an INI file which needs to be submitted as the single program parameter. The application opens a single window. Calculation starts immediately. Since waiting for the calculation result is quite boring, I built in some interactivity. So, whenever you move the mouse, a trace starting from the current mouse position is calculated and drawn into the application window. This is a toy, but can be quite useful to get an impression of the final picture by investigating certain starting points and the influence of small variations. If you came to this article looking for the butterfly effect, use this feature and watch the trace morph!

If you want to know what the application is currently calculating, press the right mouse button in order to see the traces for the pixels currently being calculated.

Application window

Image 5: Application window with a calculation in progress.

The application makes use of multicore CPUs by spawning one calculation thread per core and setting the thread affinity to this core.

Config file format

The configuration file format is the same as for all Windows INI files. You should be familiar with INI file sections, keys, and values.

Section Key Description
[FIELD] COLS Number of columns used for the field discretization. Use this parameter to set the output image width.
HEIGHT Number of rows used for the field discretization. Use this parameter to set the output image height.
SIM_WIDTH
(optional) 		Width of the field to be simulated in length units (meter). Set to COLS if unspecified.
SIM_HEIGHT Length of the field to be simulated in length units (meter). Set to ROWS if unspecified.
WIN_WIDTH
(optional) 		Width of the output window in pixels. Use this parameter to adjust the output window width when calculating images larger than the screen. Set to COLS if unspecified.
WIN_HEIGHT
(optional) 		Height of the output window in pixel. Use this parameter to adjust the output window width when calculating images larger than the screen. Set to ROWS if unspecified.
[SIMULATION] THREADS
(optional) 		Defines the number of threads spawned for the calculation. If unset, this number equals the number of processors reported by the system. Each thread runs on a different core, resulting in increased performance on multicore systems (ca. 40% on a dual core CPU).
FRICTION Friction coefficient. Friction force is proportional to this value minus velocity of the pendulum.
PEND_HEIGHT Height of the pendulum plane over the magnet plane in pixels.
DELTA_T Size of time integration steps in time units. Smaller is more accurate; this is parameter h in the Beeman integration formula.
MIN_STEPS Number of minimum steps each trace must have before an abort condition is checked.
MAX_STEPS Number of maximum steps that a trace can have. If this number is reached, the simulation stops even if pendulum did not stop over a magnet.
ABORT_VEL If pendulum drops below this velocity, it is considered to have stopped.
COLOR_SCHEME An equation used to determine the color scaling of a pixel. Color scaling depends on the trace length and the maximum trace length. This must be a mathematical expression containing the variables len and max_len (see Image 2). Example: 1/(exp(0.000001*(len*len)))
BATCH_MODE
(optional) 		Set this to one for activating the batch mode. In batch mode, the application terminates itself after finishing a calculation. This can be used for creating animations using shell scripting techniques.
[SOURCE ...]
... standing for an index 		TYPE Type of the source. Use either INV_SQR or LINEAR. Force caused by a source is either proportional to the distance of the pendulum or to the inverse square of the distance. The first is a source analog to Hooke's law, the second is a source analog to Coloumbs law (although it is not Coulombs law since I'm dealing with magnets, not electrical charges).
COLOR Color of the source.
RAD Source positions are provided in polar coordinates. This is the distance of the source from the simulation center.
THETA Source positions are provided in polar coordinates. This is the angle of the source from the simulation center.
MULT Force strength multiplicator.
SIZE Size of the source. Abort conditions are only checked if the pendulum is closer than that many pixels to a source.

Examples of configuration files can be found in the data subdirectory. If you want to create your own config files, please start by modifying one of the existing files.

Related links

The following links will provide you with additional information related to the magnetic pendulum fractal or fractals in general.

Gallery

Finally, I'd like to present some images calculated with this application. The images where taken using different parameter sets with regards to magnet strength and number, as well as mount point position.

Gallery
Image 6: Gallery showing images based on different parameter sets.

History

Rev 1.00: 30/10/2006 - Initial release

Feature list:

Rev 1.01: 19/11/2006

Feature list:

Creative Commons License
Author: Ingo Berg; Except where otherwise noted, content on this site is licensed under a Creative Commons Attribution 3.0 License.