This article is based on an article named "Experimente zum Chaos" (translated.: "Experiments on Chaos") [1] from the german magazine Spektrum der Wissenschaft. The original article introduces, among other things a simulation showing the chaotic motion of a pendulum under the influence of gravity and three magnets. This article will provide source code and explanations of the science behind simulating a magnetic pendulum.

The simulation is based on computing the motion of a magnetic 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 finally came to a rest.

When running the simulation you will see that a minor variation in the starting point of the pendulum will result in a vastly different pendulum motion after some time. Scientifically spoken one can say that the result of the simulation is highly sensitive to its starting condiftion. The effect that small variations in the initial conditions of a dynamical system will result in large, seemingly unpredictable variations of the result is investigated by a special branch of mathematics, the so called "Chaos Theory".

In 1963, Edward Lorenz an american mathematician and meteorologist discovered a similar behavior when studying the influence of initial conditions on simulations used for whether predictions. He wrote:

*
One meteorologist remarked that if the theory were correct, one flap of a sea gull's wings would be enough to
alter the course of the weather forever. The controversy has not yet been settled, but the most recent evidence
seems to favor the sea gulls. [4]
*

In later speeches the "seagull" was transformed into a more poetic "butterfly" and the term "Butterfly effect" was born to describe the phenomenon.

In the next section I will describe the underlying physical model. In the end you will see that small changes in the initial conditions and the environment of the simulation can result in large variations of the result. This will ultimately lead to a certain level of unpredictability in the simulated pendulum endpoint.

The model assumes we are dealing with a magnetic pendulum that is attracted by three magnets. Each of the magnets is assigned a distinct color (red, green or blue). The magnets are located on the edges of a equilateral triangle with the pendulum mount point being at the center. The magnetic force of is strong enough to attract the pendulum in a way that it will not come to rest in the center position. The picture below shows a rendering of the model setup (Left).

Low Friction | High Friction |

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 map showing a pattern composed of red, green, and blue pixels. The applet on the right side of the above image can be used to show samle traces. A slider can be used to set up different friction values.

For obtaining better looking results it is advisable to vary the brightness of each pixel. This can be done by taking the length of the pendulum trace into account. The diagram below shows three different color functions. They take the length of the pendulum trace as their sole parameter.

As you can see, the result is quite interesting. Starting points resulting in longer traces are shown in darker colors thus adding additional complexity to the image. The sample application provided with this article application allows you to define custom expressions 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}with:

- \(\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 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; }

You might also like: