The magnetic pendulum fractal

• Download demo project - 173 Kb
• Download source - 237 Kb

Table of contents

1. Introduction
2. The Butterfly effect
   • Simulation overview
   • Simplifications
   • Governing equations
   • Code for tracing the pendulum
3. Using the application
   • User interface
   • Configuration file format
4. Gallery
5. References
6. History

Introduction

A few years ago, I got 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 others 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 want a high resolution image of this fractal for my wall, I had to recode it. (Stupid reason isn't it? So what!)

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:

• An interest in chaos theory.
• A fast processor.
• Hardware OpenGL support. (No it's not 3D, I recycled 2D routines from my wator screen saver.)
• Quite some time.

What you get:

• A toy for creating pretty pictures and for playing
• Configure models with any number of sources using INI files
• You can abort and continue calculations

OK, since this is CodeProject, 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.

For more details, I refer you to the original: Butterfly effect at Wikipedia.

Simulation overview

This article describes 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.

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:

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. Those formulas will be interpreted using muParser, one of my other projects. (Finally, I found an application for my own project!)

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.

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.

Governing Equations

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.

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:

Code for tracing the pendulum

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

The graphical user interface

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.

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. This approach uses the Windows API rather than OpenMP which would have been easier but requires Visual Studio 2005.

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	Height 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.

• Butterfly effect at Wikipedia - Read more about the basic idea behind the Butterfly effect.

• Magnetic Pendulum Fractal Simulation On Parallel Computers - Don't want to wait for hours? This is your chance!

• Paul Nylander's Science and Art page - Mathematics can be quite beautiful, see why! (Have a look at the PovRay animation of the pendulum.)

• Urijah Kaplan page dedicated to the magnetic pendulum fractal - very informative with many case studies.

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.

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

History

Rev 1.00: 30/10/2006 - Initial release

Feature list:

• Support for a single calculation thread.
• Abort and continue.

Rev 1.01: 19/11/2006

Feature list:

• Support for distributed calculation on multiple core CPUs added.
• Creating images larger than the screen is now possible.
• Setting output window size and simulation added.
• The binary archive is now significantly faster due to using VS2005.
• Batch mode added.