Creating a pendulum simulation in C#

Question

I'm working on a game that involves swinging around a pole, and I want to simulate the physics, not just hard code the rotation. I figure a pendulum is a fairly decent model to start with, but I coded a fairly basic simulation and when I run it it slowly starts to drift. Here's a rough summary of the code I have:

    //let g be the force due to gravity
    //let o be a vector from the rotation point to the object’s current position (offset)
    //let m be the normalized cross product of o and the object’s right hand direction (aka the direction of movement)
//get the component vector of gravity in the movement direction

$$\vec h = (\vec g \cdot \vec m)\vec m$$

    //let v be the object’s current velocity and w be the component in the direction of movement

$$w = \vec v \cdot\vec m$$

   //then let c be the centripetal force

$$\vec c = -w^2 \frac{\vec o }{\lVert \vec o \rVert}$$

    //add the acceleration to the current velocity (multiplied by a small time interval dt)

$$\vec v = \vec v + (\vec c+\vec h)\ dt$$

    //change position based on velocity

$$\vec x = \vec x+ \vec v\ dt$$

If I run it for a little while, the pendulum starts slowing and drifts downwards. I can fix the drift by just hardcoding it to stay within a certain range, but I'm not sure why the pendulum is slowing.

Can anyone let me know if I'm doing something wrong, or if there's a better way to do this?

Answer 1

Numerical integration of the equations of motion of a pendulum faces two separate but interconnected problems.

The first problem is common to every Newtonian dynamical system. It is connected to how the unavoidable inaccuracy introduced by the numerical integration modifies the qualitative and quantitative features of the exact solution. In particular, it is of the utmost importance to understand the effect of the algorithmic error on the conserved quantities. In general, if an algorithm based on a fixed time step $\Delta t$ has a global error proportional to $\Delta t^n$, the energy can be conserved only at the same order. However, this is not the whole story because, depending on the algorithm, errors on a periodic motion may or may not compensate.

The simplest possible algorithm, the explicit Euler one, advances from time $t$ to time $t+\Delta t$ according to $$ \begin{align} x_{n+1}&=x_{n}+v_{n}\Delta t\\ v_{n+1}&=v_{n}+a_{n}\Delta t, \end{align} $$ but it is known to be unsatisfactory for serious numerical work (in the formulas, $a_n$ is the acceleration evaluated from the position at the time $t$). It is globally a first-order algorithm, but the conservation of energy is very poor, and even with a very small time step, results show a drift in the energy.

A much better algorithm is the Euler-Cromer algorithm, which is what has been used here. It is summarized by the evolution steps $$ \begin{align} v_{n+1}&=v_{n}+a_{n}\Delta t\\ x_{n+1}&=x_{n}+v_{n+1}\Delta t \end{align}. $$ It is still first-order but is much more stable than the Euler algorithm. The global deviations from energy conservation remain $O(\Delta t)$, but they oscillate, and there is no systematic drift of the energy.

Interestingly, without significant additional work, in the case of forces depending only on the position, one could use a definitely better algorithm, the Störmer-Verlet: $$ \begin{align} x_{n+1}&=x_{n}+v_{n}\Delta t+\frac12 a_n \Delta t^2\\ v_{n+1}&=v_{n}+\frac12 (a_{n}+a_{n+1})\Delta t, \end{align} $$ that is a second-order and symplectic algorithm.

However, this first part on the numerical integration algorithms is only half of the story when dealing with the case of a pendulum. The second half, in a way much more important, has to do with the way the constraint of a fixed length of the pendulum is incorporated in the description of the motion.

The simplest way and I would recommend it, is by recasting the description and the equation of motion in terms of a single angular coordinate. It requires to work in terms of angular velocity, angular acceleration, and torque, expressed as a function of the angle, but the final equation is pretty simple: $$ L \ddot \theta = -g \sin(\theta), $$ and can be integrated in the same way as the equations of motion written in cartesian coordinates. The main advantage is that the constraint of a fixed length of the pendulum is built-in. Obtaining the cartesian coordinates from the length and the angle is trivial.

The alternative is to work with two degrees of freedom (for instance, the two cartesian coordinates of the center of mass of the pendulum. In this case, one has to consider the presence of the reaction force due to the constraint of a motion on a circle. From the analytical side, it is not difficult to obtain the expression for such a constraining force. However, on the side of numerical integration of the equations of motion, a constraint of this kind poses a fundamental problem. The numerical evolution of the system is accurate only at the order $\Delta t^n$ ($n$ depending on the algorithm). This implies that also the constraint will be fulfilled only at the same order. And this may be a big problem. For instance, in the pendulum case, if the length varies systematically, drifting from the starting value, the period of the pendulum will be affected.

Therefore, in the case of a constrained dynamics (with this kind of so-called holonomic constraints), the numerical algorithms must be modified. The known solution of the problem (J.P. Ryckaert, G. Ciccotti, H.J.C. Berendsen. J. Comput. Phys., 23 (1977), p. 327) requires to add to the ($n+1$)-step positions and velocities an additional modification enforcing the constraint. Such modification depends on the algorithm and usually requires the numerical solution of a system of equations at every step.

It would be possible to say more on the algorithms for constrained evolution, but I think this information could be enough for the original question. In particular, as already said, I would strongly advise using the angular coordinate description.

Answer 2

For computational integration of ordinary differential equations you need a suitable time integrator. Without it, the result is likely to quickly violate energy conservation, constraints ("drifting away"), and show artificial damping. I would recommend at least that you read the chapter "Integration of Ordinary Differential Equations" from numerical recipes.

What you do by writing

v=v+(c+h)∗dt

position=position+v∗dt

is Euler integration, the most elementary of techniques to integrate ODE's.

Answer 3

Simulating physics with numerical computations is very different from working out motion with paper and pencil. First, you are dealing with computer numeric types which have finite precision where every calculation has rounding errors that can accumulate over time. Second, and more importantly, doing simulations in discrete time steps introduces errors for the same reason that numerical integration results in errors: you are replacing the actual forces and motion involved with simpler functions that can be calculated quickly and easily.

In your case, you simplify the physics by assuming that the acceleration of the pendulum bob is constant over the entire time step. This is not true of the pendulum, because the force is constantly changing as it rotates. The actual centripetal force is changing from the beginning of the time step to the end. The two expressions that update the velocity and position do not take this into account.

As an explicit example, start the pendulum at a horizontal position with zero velocity. Then, $$h=g$$ $$\omega=0$$ $$c=0$$ $$v = (0,0) + (0,g)*dt = (0,g*dt)$$ $$p = (R,0) + (0,g*dt)*dt = (R,g*dt)$$ Where $p$ is the position, $R$ is the radius of the pendulum, and every pair $(x,y)$ is a two-dimensional vector. The important thing to notice is that the pendulum bob is now farther away from the pivot. Instead of moving on a circular arc, the initial motion was straight down. This is where your downward drift comes from. The way the new velocity is computed always leaves the pendulum bob outside the radius of the current swing of the pendulum.

For every bit of physics that you want to keep consistent with reality, you need an approximation scheme that enforces it, whether explicitly in the code you write or implicitly based on the algorithm you choose. Here's an explicit method that adjusts the new position and velocity to keep the physics right.

0. At the beginning of the simulation, calculate $E = \frac{1}{2}mv^2 + mgy$. This is the total energy of your pendulum that should remain constant through the whole simulation. Then, for each time step:
1. Compute the new position and velocity as you have.
2. Scale the position vector to the actual distance from the pivot: $p = R*\textrm{unit}(p)$, where $\textrm{unit}()$ is a function that returns the unit vector in the direction of its argument (I'm assuming the pivot is at the origin).
3. Adjust the velocity so that energy is conserved. Find $v$ such that $E = \frac{1}{2}mv^2 + mgy$ continues to hold at the new $y$ coordinate of the position vector.
4. Remove any component of the new velocity that is parallel to the position vector.
5. Scale the velocity vector so that the magnitude of the vector is equal to the $v$ computed in step 3.

In short, use your initial method to approximate the new position. Then, fix up the position so the swing stay circular. Finally, adjust the velocity so the swing is circular and energy is conserved.

Now, this method has problems of its own. It probably underestimates the distance that the pendulum travels due to the scaling that puts the pendulum bob at the correct distance from the pivot. The effect of this is probably the period of the pendulum being too long for a given length. Instead of simply retracting the pendulum straight towards the pivot, you could perhaps calculate how far the pendulum bob will travel, and then put the bob at a distance over an arc. You'll have to experiment in your game to see if what results is accurate enough.

Answer 4

Most of the answers here try to deal with integrating the position/velocity of the object and subject it to the pivoting constraint. Additionally, the treatment of an object as a point mass maybe lead to unrealistic results.

My suggestion is to define the position of the object in terms of a single degree of freedom $q$ (like the rotation angle) and form the equations of motion in terms of the DOF variable(s) and their derivatives. This is called kinematics and it should be part of any simulation as it establishes what variables are needed to fully describe the configuration of the system.

For example, body (1) below is pivoting about point A and the location of the center of mass is given as a function of the angle $\theta$.

$$ \pmatrix{ x_1 = x_A + \ell \sin \theta \\ y_1 = y_1 - \ell \cos \theta } $$

Where $\ell$ is the distance between the center of mass and A, and the angle $\theta$ obeys some convention in terms of sense of rotation and where 0 is.

It follows then that the velocity of the center of mass is

$$ \require{cancel} \pmatrix{ \dot{x}_1 = \cancel{\dot{x}_A} + \dot{\theta} \ell \cos \theta \\ \dot{y}_1 =\cancel{ \dot{y}_A} + \dot{\theta} \ell \sin \theta } $$

And the acceleration

$$ \pmatrix{ \ddot{x}_1 = \ddot{\theta} \ell \cos \theta - \dot{\theta}^2 \ell\sin \theta \\ \ddot{y}_1 = \ddot{\theta} \ell \sin \theta + \dot{\theta}^2 \ell \cos \theta } $$

At this point, you apply the laws of dynamics relating forces and torques to the acceleration of the center of mass

$$ \pmatrix{ A_x + F_x = m \ddot{x}_1 \\ A_y + F_y = m \ddot{y}_1 \\ - \ell ( A_x \cos \theta + A_y \sin \theta ) &= I \ddot{\theta} }$$

where $A_x$ and $A_y$ are the pivot forces, and $F_x$ and $F_y$ the applied forces to the body (such as gravity). Also, $I$ is the mass moment of inertia of the body at the center of mass, about the pivot axis.

The solution is

$$ \ddot{\theta} = \frac{ \ell ( F_x \cos \theta + F_y \sin \theta)}{I + m \ell^2} $$

You simulation would integrate $\ddot{\theta}$ to update velocity and angles

$$ \pmatrix{t \\ \theta \\ \omega } \leftarrow \pmatrix{t \\ \theta \\ \omega } + \Delta t \pmatrix{ 1 \\ \omega \\ \ddot{\theta} } $$

or more elaborate integration schemes such as the midpoint or 2nd order RK methods.

Answer 5

thanks for all the help! In case anyone's curious, I tried testing a few different methods and you can see the results here:

The one with the red line is my original, which as you see slows and drifts.

The one with the green line is based on Mark's suggestion of trying to maintain a constant distance and energy. As you can see, it maintains distance and energy, but for some reason the direction of velocity starts to drift (I may keep working on it to see if there's a better way to maintain velocity).

Finally the one with the white line is based on Giorgio's suggestion of using angular coordinates. As you can see it works the best, though it may require more work to get it to work nicely with other physics systems in the game.