Advertisement

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.

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.

- An interest in chaos theory.
- For understanding the code i recommend a basic understanding of C++ and MFC.
- Hardware OpenGL support is required for the rendering.

- A software for creating pretty high resolution images
- Configure models with any number of magnets using INI files
- You can abort and continue calculations and make use of multiple processors for speeding up the computations.

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]:

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:

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.

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:

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.

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.

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.

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}wobei:

- \(\vec{a}_t\) - Total acceleration of the pendulum
- \(\vec{a}_g\) - Acceleration caused by the pullback of the pendulum
- \(\vec{a}_m\) - Accelleration caused by the magnet m
- \(\vec{a}_f\) - Decelleration caused by friction
- \(k_f\) - Coefficient of friction
- \(k_g\) - Strength of the pendulum pullback
- \(k_m\) - Strength of the the magnet
- \(\vec{r}\) - Position vector of 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; }

Advertisement

You might also like: