consider modeling a double pendulum or its generalization to more segments with a numerical ordinary differential equation solver. (the generalization does not seem to have a name: the term "compound pendulum" refers to something else. perhaps call the generalization "a multiple pendulum", not to be confused with "multiple pendulums" which is just a collection of pendulums / pendula.)
motion of a frictionless (undamped) double pendulum makes intuitive sense but looks aesthetically bizarre.
Verlet method, Stoermer's rule, and symplectic methods initially look promising because they are numerically stable. Runge-Kutta is not numerically stable.
numerically stable does not mean that energy is conserved. it only means energy does not head to infinity. in practice, for Verlet, energy oscillates around the starting value. however, i suspect that, because of roundoff error, even for numerically stable methods, energy can eventually differ from starting energy by an arbitrary large amount if you let the simulation go for long enough, behaving like a random walk.
because effects like cracking the whip can occur, if we want accuracy, we need to use an adaptive stepsize algorithm. but maybe we don't need accuracy if it only has to look good enough for entertainment (computer graphics).
adaptive stepsize Verlet methods are few and far between. the only one i could find is in Numerical Recipes c16-5, "Second-Order Conservative Equations", which plugs Stoermer's rule into the Bulirsch-Stoer method. adaptive stepsize for Bulirsch-Stoer is described in the previous section, c16-4.
Peter Young's notes, "The leapfrog method and other "symplectic" algorithms for integrating Newton's laws of motion" http://physics.ucsc.edu/~peter/242/leapfrog.pdf, curiously state "I am not aware of adaptive stepsize control for symplectic algorithms" one sentence after mentioning Numerical Recipes.
the notes also disrecommend regular Verlet because it is more susceptible to roundoff error than velocity Verlet. Numerical Recipes's Stoermer's rule mentioned above gives a method (from Henrici) to reduce that roundoff error.
despite these criticisms, the notes are pretty good, giving the best explanation i could find for Verlet's numerical stability (computing the Jacobian), and explicitly saying energy is not conserved, but listing other things which are conserved. it also gives higher order symplectic algorithms.
incidentally, we also need adaptive stepsize for N-body gravity simulations because close passes can occur. this and crack-the-whip mentioned earlier seemingly eliminate the naive (constant-stepsize) Verlet method from the domains it is most famously designed for.
two constraints: pendulum segments are fixed length (rigid). conservation of energy.
keeping pendulum segment lengths constant requires some sort of tweak no matter what algorithm we choose. unless we really don't care about accuracy: we could let the segment lengths drift over time. we won't explore this further.
let segment positions be in Cartesian coordinates. the force of gravity acts on the center of mass, while the segments act on each other through pivot points. we briefly thought about using polar coordinates, defining a local coordinate origin for each segment at its pivot point, and doing a coordinate transform for each segment going down the line, because that seems complicated. however, others have done it, via Lagrangians: https://jakevdp.github.io/blog/2017/03/08/triple-pendulum-chaos/
probably bad idea: instead of being rigid, model segments as very stiff springs. unfortunately, the stiffer the spring, the smaller the step-size you need for accuracy. maybe it's possible to do a Richardson extrapolation-like trick to extrapolate to infinite stiffness each time step. once a spring starts vibrating, it's going to keep vibrating forever because we have no damping. if we were to add damping, then energy would no longer be conserved.
not so bad idea: let the segments be not-very-stiff springs. this completely changes the model, but it may be acceptable, even desirable, for entertainment: jiggly pendulum. it's curious how making the model more complicated might cause solving it to become more straightforward.
tweak for rigidity: take one time step with the ODE solver. move the endpoint of the first segment horizontally until its length is correct. then move the endpoint of the next segment horizontally to meet the rigidity constraint, and so on down the line. nothing changes height, so gravitational potential energy remains unchanged. (this assumes segments have constant linear density.) there are generally two horizontal positions which satisfy rigidity: pick the closer one. there may be some tricky details. even more tricky if 3D.
tweak velocities so that energy remains conserved each time step. one cannot use regular Verlet for this because we need velocities to calculate kinetic energy.
more complicated: adjust velocities and positions simultaneously, a constrained minimization problem, probably non-linear. adjust both velocities and positions as little as possible (by some measure) to meet the constraints. should we try to minimize adjustments to velocity or to kinetic energy? they differ by a square root.
given that we are tweaking for rigidity and conservation of energy every time step, stability of Verlet no longer provides much of a benefit.
can one do similar tweaks to conserve energy in a N-body gravity simulation? also conserve momentum.
segments can be decorated to be any shape. center of mass of a segment can be anywhere: compound pendulum. interesting might be tracking a point on the opposite side of the pivot point from the center of mass.
No comments :
Post a Comment