There's a terrific test of the accuracy of a pendulum simulation, one that many workers overlook. The test is based on the fact that, for a swinging pendulum, gravitational potential energy (GPE) and kinetic energy (KE) must sum to a constant everywhere in the arc. The reason is that energy is conserved. So if a simulator shows a deviation from this constant sum, there's something wrong with the simulator and/or the underlying mathematics.

This thinking also applies to orbital simulations -- in an elliptical orbit, an object's radius from the central mass, and its velocity, are constantly changing, but (barring friction) its energy, the sum of GPE and KE, must be a constant.

> So if a simulator shows a deviation from this constant sum, there's something wrong with the simulator and/or the underlying mathematics.

Not necessarily so -- many (most) discrete integrators do not exactly conserve the Hamiltonian. In fact, large classes of geometric integrators yield noticeable (but bounded) amounts of 'noise' on top of the true Hamiltonian.

>> So if a simulator shows a deviation from this constant sum, there's something wrong with the simulator and/or the underlying mathematics.

> Not necessarily so

But the "constant sum" I refer to is the sum of gravitational potential energy and kinetic energy in the system. If that were not true, if the sum were not a constant (and remembering this is an ideal model without friction), then energy would not be conserved. So, yes, necessarily so.

Is it possible to apply that test for a multi-body simulation? I guess you could look at the combined GPE + KE of the entire system, but not individual bodies?

You can apply it for any system with only conservative forces. And gravitation is conservative.

https://en.wikipedia.org/wiki/Conservative_force

This would need way more explanation to make sense. The part I like is that the equation pyramid looks like a Christmas tree ;)

The principle of least action states that the path that a system takes is the one that minimizes it's 'action'. Famously, this says that light travelling through a lens will go through in the fastest path (taking into account the fact that it is slower moving through the lens than through the air). The Euler-Lagrange equations give an approach to solving 'any' classical system using this principle. In order to do so, you first parametrize the system by specifying some variables that describe the state of the system. Then you describe the energy of the system in terms of its parameters. Then you plug that into Euler-Lagrange and you can solve for the motion of the system based. (The 'action' in this sense is determined by the energy along a path that the system might take.)

(This is in contrast to Physics 101 approach of a free body diagram where you add up a bunch of forces on an object. In particular the 'constraints' that the system faces (in this example, it's that the pendulum is confined to moving along a circle) are implicitly included in the parametrization of the system and so there's no need to ever consider directly the force of the bar of the pendulum that holds the weight in place.)

The first example on the page parametrizes the pendulum in the 'obvious' way, by noting that you know the state of a pendulum if you know what angle the weight forms relative to the vertical (plus also the rate at which that angle is changing). This is nice, but the solution involves trigonometric functions.

The second one parametrizes the same system in a different way. Here the mu parameter is, roughly speaking, the height of the weight above the minimum (but not really). It still specifies where the pendulum is (the equation M(mu) = ... tells you the (x,y) coordinates of the pendulum's weight for a particular mu), so we can still calculate the energy of the pendulum in terms of mu. Energy has two components, kinetic and potential. Potential is the just height from 0 times mass times gravitational constant g. The kinetic is mv^2/2 and the speed v can be determined by using the parametrization above M(mu)= .... The end result is a different means of expressing the same system, but interesting in a manner that does not involve trigonometric functions despite being 'so obviously' trigonometric in nature.

One might hope that eliminating a trig function would speed up simulation. I think it's clear that one would need a more demanding simulation to make the difference noticeable.

I think this is a fun little project and it was nice to see it introduce me to a new concept (rational trig).

If you are using floating point numbers, you are using rationals. You cannot represent an irrational number using a floating point variable. There is no getting around the fact that the real number system is an abstraction, and floating points model a very small segment of the reals (the rationals are of measure zero, and hence so are all the reals representable as floating points). There are only so many bits with which to do computations, according to the IEEE 754 scheme.

Edit: addendum

How exactly does one do ordinary calculus with only rational numbers?

The use of the rational numbers I was referring to was more in the equations than in the running code.

In fact, part of the appeal I find in rational trigonometry is that it should be easier to adapt for computing, since as you pointed out, computers can not represent irrational numbers. Yet trigonometric functions usually return irrational numbers, and a computer can't do better than approximate them.

On the contrary, it is possible to do exact arithmetic on rational numbers. Granted, you can't do that with floating points, though.

PS. To be clear: although the numbers actually used in the code are floating point numbers, the equations have a purely algebraic form (no sine, cosine or sqrt), so it should at least in principle be possible to use a Rational number type instead, and the computation would then be exact. The only approximation would be due to the Runge-Kutta algorithm.

On the contrary - computers certainly can represent irrational numbers - but not with floating points.

Computers, like humans, can only represent a vanishingly tiny proportion of “real numbers”, either by giving specific ones a name, or some kind of symbolic formula by which to compute with them. To do practical computations we need to break everything down into basic arithmetic operations, which ultimately means we have a concrete multiplication/addition(/other operation) table encoded somewhere, e.g. in silicon.

One that I like for doing exact arithmetic on numbers represented by black-box functions is a continued fraction representation, see http://perl.plover.com/yak/cftalk/INFO/gosper.txt

(Note: floating point is successful because trying to use irrational numbers in computations is most of the time a huge waste of computational resources compared to approximating with a binary fraction.)

Here’s an interesting precise representation of algebraic numbers as roots of integer polynomials: http://twistedoakstudios.com/blog/Post6871_impractical-exper...

That's pretty much how a typical CAS represents algebraic numbers, I don't think there's much impractical about that idea.

No, but algebraic numbers are only a tiny portion of all irrational numbers and importantly misses key ones like e and pi, so the original point still stands.

The rationals may be measure 0 but they do still approach the irrationals to any desired degree.

The rational parametrization looks like the Cayley transform, which IIRC does not cover the whole circle (you can't reach (-1, 0) with a finite mu). Does this happen here as well?

Indeed, (-1, 0) can not be reached. I'm thinking of ways to improve it that regard.

This is a stereographic projection representation of an angle. The number should be seen as a point on the real projective line, which includes a point at infinity.

If you’re using floating point, you can use the number ∞ for this. You’ll need some special cases in the code for it.

Cf. https://en.wikipedia.org/wiki/Stereographic_projection

https://en.wikipedia.org/wiki/Tangent_half-angle_substitutio...

https://en.wikipedia.org/wiki/Real_projective_line

Very neat! But possibly afflicted somewhat with rounding errors; I spent a few minutes reading about Lagrangian mechanics, and came back to find the pendulum whirling around its origin point like a thrown bolas. It looks like you've put the proper amount of effort into scope hygiene and the debugger doesn't seem to want to work, so I can't give the values of any of your code's variables, but I'm using Firefox 45 ESR, if that helps at all.

I've also noticed that the simulation goes berserk, but only if you leave the page (I mean if you go to an other window or tab). I don't think it's due to Runge-Kutta but rather to the way HTML5/javascript behaves when the page is running in the background.

I have no idea what could be done to avoid it.

This happened to me, too. My guess is that the simulation stops running when not open, but when it next runs it computes using a time delta calculated by the actual time since the tab was seen. So the integration adds a ridiculous amount of acceleration to the velocity, at which point the whole thing is doomed. This problem: http://gafferongames.com/game-physics/fix-your-timestep/

Yep. This is the code that's causing it[1]:

  var now = Date.now()/1000;
  var dt = now - then;

Adding this line fixes it:

  dt = dt > 1 ? 1 : dt;

[1] https://github.com/grondilu/grondilu.github.io/blob/47d41d/p...

I bet you're right! Or close, anyway. From [1]:

> The callback rate may be reduced to a lower rate when running in background tabs or in hidden <iframe>s in order to improve performance and battery life.

[1] https://developer.mozilla.org/it/docs/Web/API/window/request...

That looks like a sensible explanation, thanks. I guess I can put a ceiling on dt as a quick work-around.

Cool stuff! It's always nice when you can formulate your problem in a way that avoids transcendental numbers and functions.

Your page contains a typo, however: “The Euler-Lagrance equation then gives...”

Please change the name then to "rational trigonometry", not rational numbers. I was thinking that the only numbers you were using involved the ratio of integers, but I did not see this explicitly here. That or I don't read it correctly.

There's a specific definition of the set of rational numbers which doesn't quite mean what I think you meant when you include it in your title. I mean, it's your code, you decide what you call it but it will most likely confuse others.

Ok, we s/numbers/trigonometry/'d the title.