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.

You must Sign In to use this message board.

FAQ    Noise Tolerance Very HighHighMediumLowVery Low   Layout Expand Posts & RepliesThread ViewNo JavascriptNo JS + Preview   Per page 102550

Msgs 1 to 14 of 14 (Total in Forum: 14) (Refresh) FirstPrevNext

Great going! Excellent effort Yogesh P. Dhakad 7:31 17 Feb '07   Keep it up man. You get my five points. Regards, Yogesh Dhakad Sign In·View Thread·PermaLink Awesome! Tydren 19:13 22 Nov '06   Dude, this is truly some awesome stuff! Excellent work, you got my 5! Sign In·View Thread·PermaLink Article Updated to V1.01 iberg 11:20 22 Nov '06   The article has been updated. New is the support for multiple core CPU's and improvements in the image save routines that now instead of taking screenshots create the bmp file in memory. Hence images can be larger than the screen now and windows do not have to be visible when the image is saved. Sign In·View Thread·PermaLink It is NOT a magnetic pendulum KevinMcCann 23:55 8 Nov '06   Coulomb's law is the electrostatic (not magnetic) interaction between two charged bodies. The magnetic interaction would be approximately that between two magnetic dipoles and is more complicated and behaves as 1/r^3. People use the electrostatic interaction because it is easier to program, i.e. most of the sites that purport to do the magnetic pendulum are done by people who either don't know or are too lazy to do it right! The demo is nice, but the title is incorrect. Why not change it to the electrostatic pendulum and at least get that part right? Kevin Sign In·View Thread·PermaLink 1.00/5 (1 vote) Re: It is NOT a magnetic pendulum iberg 9:44 9 Nov '06   There is no need to change anything in the article since you do not seem to have read it. Had you read the article you would have found out rather quickly that the coulomb law is only used to illustrate the type of force, you probably would have become suspicious when reading about the law of gravity too in this context, you would have read that there is mention of magnets actually being dipoles and the force for magnets being proportional to 1/r^3. And you would have found my question asking for explanation why this fact is ignored often normally. But maybe you read the article and this is your answer. So lets consider it: Your statement it is more complicated to calculate dipole forces is incorrect. Assuming a force proportional to 1/r^3 requires changing on byte in the code assuming the pendulum is in the far field of the dipole which is reasonable since the magnets are flat. The dipole field is just the result of a superposition of two monopol fields, one positive the other negative. This results automatically in a field proportional to 1/r³ assuming you are in the far field of the dipole. Even when assuming you are not in the far field the calculation does not get more complicated. The only thing you have to do is to position three additional sources with negative coefficients directly below the three sources used in this simulation. Two monopol sources with opposite coefficents placed on top of each other yield what? Right a Dipole! It would probably slow down the calculation due to a doubling of involved sources it does not increase the complexity of the problem. There is one case that could potentially increase the complexity namely assuming the pendulum itself is a magnet but as far as this simulation is concerned the pendulum is made up of iron. Which means it is not a dipole in itself. When dealing with dipoles it is reasonable to use forces proportional to 1/r² when you are close to one of the poles because the field of each pole is proportional to 1/r² and the other pole is too far in order to have significant influence. For this simulation this would mean the magnets are actually bar magnets and they are longer than the distance from the magnet plane to the pendulum plane. Call the forces electrostatic, call it gravitational call it a hypothetic magnetic monopol or what ever you like. It does not matter. It is all equally correct or equally wrong. As far as the model is concerned. This software simulates pendulum consisting of a infinitively small iron mass hanging on a massless stick of infinite length over a plane being attracted by three bar magnets extending from the pendulum plane downwards to infinity. Sign In·View Thread·PermaLink 2.00/5 (1 vote) Re: It is NOT a magnetic pendulum KevinMcCann 23:15 9 Nov '06   Actually the interaction between a magnet and a small iron mass is a dipole-induced dipole interaction that falls off as 1/r^5. Even the dipole-dipole interaction is more complicated than you make out, since the relative direction of the two dipoles with the vector from one to the other must be taken into account. To a first approximation, two magnets side by side with their dipoles pointing up, say, do not interact; whereas, if one is above the other, there is a strong interaction. Anyway, all that aside, it is a very nice demo of C# and scientific programming. Examples of such programming are not so easy to find in C#. Any language that requires an explicit function call to square something, pow(x,2), instead of x^2 or x**2, was not thinking about those of us who do science. Bill Gates et al had a chance to fix this oversight in C#, but did not. However, I do enjoy C# a lot more than C++ or (ugh) C. Sign In·View Thread·PermaLink Re: It is NOT a magnetic pendulum iberg 9:38 10 Nov '06   The iron pendulum will move throgh a magnetic field so eddy currents are inducted they themself create a magnetic field that interact with my field and the assumtion of a static field is gone. If this is what you mean I totally agree. Nevertheless the orientation of the dipoles is taken into account automatically by the position of the monpol sources that form them. No vector algebra involved. Just plain field theory since any potential field and we are talking about a potential field is formed by a superposition of multiple (monopol) sources. Vector algebra starts when taking the gradient of the field for the force calculation. If you like C# you may like Java too and even if not i could recommend you this page: http://www.falstad.com/mathphysics.html It is fully loaded with excellent interactive physical examples. Just awesone... Sign In·View Thread·PermaLink Awesome! jerry_pendergraft 9:36 7 Nov '06   Those images are killer man. I'm looking forward to creating some myself. Sign In·View Thread·PermaLink Texture Mapped Sphere John Stewien 1:51 7 Nov '06   Nice Demo, I assume the ultimate result would be a texture map for a sphere showing all the possible starting points for the pendulum. That Beeman's algorithm looks very useful. Update: Oops, just noticed you mentioned the sphere texture map thing. There's lots of interesting things to try out though with these algorithm. I'll have to play around a bit. Sign In·View Thread·PermaLink Re: Texture Mapped Sphere iberg 8:48 7 Nov '06   The mapping to the sphere needs a different implementation since it must use spherical coordinates. It's not that complicated once you have a vector class and the basic vector algebra implemented. Extending this to 3D requires just adding another dimension and Mapping the accelleration vectors to a sphere. I wonder why no one did this before. I don't think I'm going to do it either, well not immediately... Sign In·View Thread·PermaLink Excellent demo Duncan12345 22:40 6 Nov '06   This is a superb demo. In descriptions of chaotic patterns, the pendulum is often described, but one so rarely sees a coded example (they tend to go for complex root of -1, which can look similar but is easier to code!). And as if that wasn't enough, you've added the interactive aspect with the pendulum's path being plotted from the mouse position. I love it! Sign In·View Thread·PermaLink Re: Excellent demo iberg 8:48 7 Nov '06   Thanks a lot. Great you like it! Sign In·View Thread·PermaLink Nice Ralph Wetzel 9:31 1 Nov '06   I didn't look @ the implementation, but the article and the pictures are really nice. Got my 5! Greetings, Ralph Sign In·View Thread·PermaLink Re: Nice iberg 22:18 2 Nov '06   Thanx, you're welcome. Sign In·View Thread·PermaLink

Last Visit: 16:52 9 Jul '09     Last Update: 16:52 9 Jul '09 1