# Creating a pendulum simulation in C#

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?

• Implementation details of computational tasks are off-topic on this site: “While computational physics is on topic, we are not a programming site. If your question is about implementing computational code - in particular, if it's about writing, compiling, debugging or optimizing code, or about a specific language or library - then it is off topic.” – G. Smith Jun 6 at 3:42
• Related: physics.stackexchange.com/questions/404790/… and other Linked questions. – Emilio Pisanty Jun 6 at 4:38
• @EmilioPisanty the problem of simulating the motion of a pendulum is not the same as the simulation of the solar system. There is an additional difficulty related to the presence of a rigid constraint. – GiorgioP Jun 6 at 4:43
• @GiorgioP I think you misunderstand me. If the OP were to actually ask about the conceptual confusion then this would be a valid post; I'm not saying the issue is in the content/topic of the question. Just like how in other check my work questions we say to actually focus on the conceptual confusion. – BioPhysicist Jun 6 at 15:21

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.

• Thank you for the thorough, well explained response. And yeah, that makes a ton of sense. I think I was trying to avoid using angles since this will be in 3D and I didn't want to figure out how to calculate the angle in 3D space, but it sounds like it might be a lot simpler than what I'm trying to do. – Adam L. Jun 6 at 17:29

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.

• The algorithm is not the explicit Euler. It is the semi-implicit Euler algorithm, also known as Euler-Cromer. It has a much better behavior than the explicit Euler integration. – GiorgioP Jun 6 at 5:03
• Interesting. I’m using Unity, which has a variable called deltaTime which is the time between frames. Are you saying that it’s too big of a dt to be accurate? – Adam L. Jun 6 at 5:11
• @AdamL. When you approximate an integral by doing pos += v * dt that's called Euler integration, and it generally doesn't work very well (unless v is constant), no matter how small you make dt. From Emilio's link: "When it comes to these systems, the most simple updating scheme (Euler method) is inappropriate because it does not conserve energy, so one needs to rely on what are called symplectic integrators". However, as G. Smith mentioned, there are other problems with your algorithm. – PM 2Ring Jun 6 at 5:26
• @Adam One popular symplectic integrator is Verlet, but I prefer synchronised Leapfrog. The good news is that they're both quite easy to code. – PM 2Ring Jun 6 at 5:48
• @AdamL. Unity does not know how you estimate your state functions inbetween frames. So it can't tell you whether the time step is "too big". It is your responsibility to provide good enough estimates. In special cases you might achieve that by knowledge about the nature of the differential equations, but generally it is much easier to use an automatic/adaptive time integration. – oliver Jun 6 at 6:00

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.

1. 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:
2. Compute the new position and velocity as you have.
3. 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).
4. 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.
5. Remove any component of the new velocity that is parallel to the position vector.
6. 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.

• Why one should scale only the velocity? Energy depends on both position and velocity. – GiorgioP Jun 6 at 6:28
• @GiorgioP The position is modified to put the pendulum at the right radius. The velocity is scaled so that the energy is correct at the corrected position. – Mark H Jun 6 at 6:31
• I am not aware of any analysis of the error introduced by the algorithm you have described. Do you have a reference? – GiorgioP Jun 6 at 6:35
• @GiorgioP I don't have a reference since I just made up this algorithm. I added a last paragraph to point out potential problems. I doubt it holds up at all for real physics use, but it might be good enough for a game. Plus, it makes the physics problems explicit. If I were to do this simulation, it would be closer to your answer with angular coordinates, building in the constraints into the equations. – Mark H Jun 6 at 6:45
• This actually makes a ton of sense and seems like a fairly elegant and simple solution. I'll have to test it to see if it actually works but I like that it incorporates velocity into it, since it's possible you'll grab the pivot point with some initial velocity, and I want that velocity to be conserved as well. – Adam L. Jun 6 at 17:26

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.

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.

• Very interesting. I agree that Giorgio's is the simplest way to go about this. When you implemented my method, did you make sure to subtract the velocity component parallel to the pendulum arm before rescaling the velocity? The velocity should always be parallel to the pendulum arm. It looks like there are several times where the velocity vector is parallel to the pendulum arm, which is not possible. – Mark H Jun 6 at 20:30
• Ah no good point, I did not. When you say remove the parallel component, is that the same thing as isolating the perpendicular component? – Adam L. Jun 6 at 21:13
• Yes. The dot product of the velocity and the pendulum arm should always be zero. This is because the pendulum bob never changes distance from the pivot. – Mark H Jun 6 at 22:02
• Well will you take a look at that? imgur.com/ruFfV5V – Adam L. Jun 6 at 22:28
• Neat! As I thought, my code results in a pendulum that swings too slowly. The position rescaling causes the movement along the arc to be less than $v*dt$. You could do something like calculating the new position based on moving $v*dt$ along a circular arc, but at that point the algorithm would be converging on Giorgio's algorithm. – Mark H Jun 6 at 22:42