js'' has a very good point - when in doubt "Energy makes it go!". Messrs Lagrange and Hamilton have a nice angle on this but this can probably wait until you decide to do a double- or n-link pendulum.

Nevertheless, energy should feature in your simulation as it gives a clean and easy way to determine precisely the type of motion you want (whether oscillatory or rotational) and it also provides a scheme for stopping numerical error accumulating as the simulation runs.

js was quite specific but more generally the principal is (at least until you introduce friction or someone kicks the pendulum) that energy is a constant of the motion, ie

T + V = E = constant [1]

T denotes the kinetic energy of the bob, and V is it''s potential, these are respectively functions of the angular displacement x say and the angular velocity xdot - for this case we have

T(xdot) = 0.5 . m . (r . xdot)^2, [2]
V(x) = m. g. h(x) [3]

h = r . sin(x) or cos(x) or r(1 - sin(or cos)(x))
(precise choice of trig function depends where you want V = 0 and where you start measuring x from).

What this all boils down to is that if you''ve decided how much energy you want in the system, you can pick any x or any xdot (that is small enough to satisfy [1]) and then just solve for the other. This gives you a way to start the system going and also a way to correct it as you go along whilst using an Euler or any other method to "estimate" x|t+dt (x at time t + dt).

Here''s one possible scheme:

At the start of the simulation (time, t = 0), pick x|0 and xdot|0 "to taste". Then compute E by [2] + [3] = T + V = E. Assuming that each time step is of length dt, where dt is the time it takes to repaint your display (plus not much for the computation) you''ve already decided xdot|0 so you can precisely compute the acceleration at t = 0, xdotdot|0 and estimate x|dt using Euler, i.e.

xdotdot|0 = g . sin(x|0) // exact
x|dt = x|0 + ((xdot|0) . dt) + (0.5 . xdotdot . dt^2) // estimate

Above, the first computation is exact upto floating point precision, the second is a "guesstimate" using Euler. Now here''s the correction at the start of the next frame,

T + V = E
=> T = E - V
=> 0.5 . m . (r . xdot)^2 = E - m. g. h(x)
=> xdot|dt = sqrt[ (2/m.r)( E - m . g . h(x|dt)) ]

so, xdot|dt = sqrt( E - V(x|dt)) // exact

We''ve now got an exact xdot|dt, corresponding to our estimated x|dt and we''re back to our t = 0 situation waiting to compute x|dt+dt for the next frame.

No errors then?

Well not quite... Any error we make between frames will eventually accumulate over time so the simulation will loose or gain time as compared to a real pendulum (a bit like a clock), but you''d need a stopwatch or an actual pendulum sitting (in a vacuum) next to the monitor to be able to tell - locally (in time) the motion looks right - you can improve on the accuracy of your "clock" by picking a more accurate numerical integration scheme (this is your "watch mechanism") but the improvement won''t be easy to see over a few minutes.

Finally, here''s the user interaction - js told you how to set the energy level to get rotational motion ie make sure that T + V > E_rot where E_rot = 2mgr. All that remains is the user interaction. You can handle this by feeding it straight into the energy constant, eg. press "plus" to increment E, press "minus" to decrement E. Then provided you''re careful and don''t let the energy go lower than the current height of the bob at the moment the user presses "minus", it will all work seemlessly because the correction step will effect a corresponding decrease in kinetic energy and the corresponding maximum height of the bob on the next full swing. #

The problem with the energy approach is that you have to make special arrangements to decide the direction of motion. "xdot|dt = sqrt(...)" only gives the correct magnitude, not necessarily the right direction... As for numerical stability, you have to watch out for negative roots and cases where the angular velocity is very close to zero. The force/acceleration approach descibed above doesn''t have these problems...

AP makes some good points here - most importantly that the energy approach is not a solution.

In all fairness though, I don''t think I was proposing the energy approach as the solution - the energy stuff is for correcting the numerical inaccuracy in the integration step which is still the driving mechanism. In fact the effects of this correction will be hardly discernable for the first few frames. What it''s for is to ensure that the pendulum doesn''t evolve away from its original "phase trajectory", ie that if you start with oscillatory you don''t gradually evolve into rotational and vice versa by iteratively accumulating error.

Having said that, AP has a point on both other counts, the xdot|dt correction method is incomplete as stated (apologies) because of the (rather crucial) sign business; and it''s good to be suspicious of square rooting things. However, the ideas are basically sound and the algorithm is easily fixed.

First the sign of xdot...
Euler gives a first guess for x|dt, (we''re all agreed on that), but Euler can also be used to estimate xdot|dt, so it can give the required sign,

i.e.

(xdot|dt)_Euler = xdot|0 + (xdotdot|0) . dt // Euler step for sign of xdot
xdot|dt = signof[(xdot|dt)_Euler] . sqrt( E - V(x|0)) // Energy approach for mag of xdot

...that''s at most one extra line of code on the previous suggestion.


Now the sqrt...
Watching out for small and/or negative roots is a very good general principle. Sqrt is a (sometimes unnecessary) computational expense and has a domain for numerical purposes of ({0} ^ (tol,large)) obviously the result is complex for negative inputs and it also performs increasingly badly as inputs approach zero (or floating point tolerance) - but that doesn''t mean that you should never use it.

In this situation sqrt(E-V) is theoretically okay - remember that V(x) is tightly coupled to E so if E-V < 0 it can only be due to a slight numerical error (or a mistake in the logic) - consequently provided the logic is okay this term can only go very slightly negative. Okay, as far as sqrt is concerned that''s like being very slightly pregnant! Which is why this is a good point. However, again there is a simple fix - only apply the sqrt when |E-V|>tol, (ie almost always) or alternatively, take E-V to be zero at these points.

...thats another line of code reading something like
if fabs(E-V) > tol. { do correction }


The only real nasty (on the "small xdot" front) is in the case of pathologically low Energy situations - there''s an analogous puzzle with bouncing balls. If you do a bouncing ball with coefficient of restitution < 1, the ball has less energy (hence reaching a lower maximum height) after each bounce, it eventually reaches a point (in time) where it''s bouncing not very high but with very high frequency and never really stops - the usual trick is to use some form of thresholding fudge to stop the ball "somewhere sensible". That''s easily doable with the pendulum - okay, it''s a fudge but just for this very special case. If you really want to study what''s happening during low energy states this becomes a totally diffrent ball-game (so to speak) - a proper treatment of the pathologically low energy case is probably beyond the scope of the original exercise(?).

quote:

---

Original post by octal
approach as the solution - the energy stuff is for correcting the numerical inaccuracy in the integration step which is still the driving mechanism. In fact the effects of this correction will be hardly discernable for the first few frames.

---

Very good point. The ultimate goal is to find a numerical method that conserves everything it should: momentum, energy, mass (which is an issue in fluid simulations). And it can be a multi-step method that uses an energy-based correction after the time step integration.

Graham Rhodes
Senior Scientist
Applied Research Associates, Inc.

Octal, I agree. If, for some reason, you would like the simulation to go on indefinitely with the exact same amplitude, you could very well use conservation of energy principles to correct for numerical errors. To avoid the problems I pointed out above, you could simply apply this correction once per cycle, when the pendulum passes its lowest point. However, in most real world applications, you would probably want to include friction, user interaction and/or other non-conservative forces into the simulation.

Wow, this is getting very involved in deed. I''m going to try and patch some kind of calculation from all the ideas here and see what results I get. Thanks for all your help!!!

quote:

---

Original post by Anonymous Poster
To avoid the problems I pointed out above, you could simply apply this correction once per cycle, when the pendulum passes its lowest point.

---

There''s no advantage to doing it only once per cycle - you want to correct as often as poss - especially if you''re thinking of scaling up to "real-world" applications.

quote:

---

Original post by Anonymous Poster
However, in most real world applications, you would probably want to include friction, user interaction and/or other non-conservative forces into the simulation.

---

Very app dependent (I suggested a way to deal with user input) but for eg, friction causes "energy dissipation" - (decrease in energy over time/distance or whatever other coordinate seems convenient) you can still track change in energy by isolating non-cons forces and differentiating/integrating, whatever''s required... eg Fdx = workdone = energychange, Fdt = impulse = momentumchange.