Why is my code for a double pendulum returning NaN?

Question

When I print my acceleration and velocity, they both start (seemingly) normal. Shortly, they start getting very big, then return -Infinity, then return NaN. I have tried my best with the math/physics aspect, but my knowledge is limited, so be gentle. Any help would be appreciated.

float ang1, ang2, vel1, vel2, acc1, acc2, l1, l2, m1, m2, g;

void setup() {
  background(255);
  size(600, 600);
  stroke(0);
  strokeWeight(3);
  g = 9.81;
  m1 = 10;
  m2 = 10;
  l1 = 100;
  l2 = 100;
  vel1 = 0;
  vel2 = 0;
  acc1 = 0;
  acc2 = 0;
  
  ang1 = random(0, TWO_PI);
  ang2 = random(0, TWO_PI);
}

void draw() {
  pushMatrix();
  background(255);
  
  translate(width/2, height/2); // move origin
  rotate(PI/2); // make 0 degrees face downward
  
  ellipse(0, 0, 5, 5); // dot at origin
  
  ellipse(l1*cos(ang1), l1*sin(ang1), 10, 10); // circle at m1
  ellipse(l2*cos(ang2) + l1*cos(ang1), l2*sin(ang2) + l1*sin(ang1), 10, 10); // circle at m2
  
  line(0, 0, l1*cos(ang1), l1*sin(ang1)); // arm 1
  line(l1*cos(ang1), l1*sin(ang1), l2*cos(ang2) + l1*cos(ang1), l2*sin(ang2) + l1*sin(ang1)); // arm 2
  
  float mu = 1 + (m1/m2);
  
  acc1 = (g*(sin(ang2)*cos(ang1-ang2)-mu*sin(ang1))-(l2*vel2*vel2+l1*vel1*vel1*cos(ang1-ang2))*sin(ang1-ang2))/(l1*(mu-cos(ang1-ang2)*cos(ang1-ang2)));
  acc2 = (mu*g*(sin(ang1)*cos(ang1-ang2)-sin(ang2))+(mu*l1*vel1*vel1+l2*vel2*vel2*cos(ang1-ang2))*sin(ang1-ang2))/(l2*(mu-cos(ang1-ang2)*cos(ang1-ang2)));
    
  vel1 += acc1;
  vel2 += acc2;
  
  ang1 += vel1;
  ang2 += vel2;

  println(acc1, acc2, vel1, vel2);
  
  popMatrix();
}

Answer 1

You haven't done anything wrong in your code, but the application of this mathematical technique is tricky.

This is a general problem with using numerical "solutions" to differential equations. Similar things happen if you try to simulate a bouncing ball:

//physics variables:
float g = 9.81;
float x = 200;
float y = 200;
float yvel = 0;
float radius = 10;

//graphing variables:
float[] yHist;
int iterator;

void setup() {
  size(800, 400);
  iterator = 0;
  yHist = new float[width];
}

void draw() {
  background(255);
  y += yvel;
  yvel += g;
  if (y + radius > height) {
    yvel = -yvel;
  }
  ellipse(x, y, radius*2, radius*2);
  //graphing:
  yHist[iterator] = height - y;
  beginShape();
  for (int i = 0; i < yHist.length; i++) {
    vertex(i,
           height - 0.1*yHist[i]
          );
  }
  endShape();
  iterator = (iterator + 1)%width;
}

If you run that code, you'll notice that the ball seems to bounce higher every single time. Obviously this does not happen in real life, nor should it happen even in ideal, lossless scenarios. So what happened here?

If you've ever used Euler's method for solving differential equations, you might see something about what's happening here. Really, what we are doing when we code simulations of differential equations, we are applying Euler's method. In the case of the bouncing ball, the real curve is concave down (except at the points when it bounces). Euler's method always gives an overestimate when the real solution is concave down. That means that every frame, the computer guesses a little bit too high. These errors add up, and the ball bounces higher and higher.

Similarly, with your pendulum, it's getting a little bit more energy almost every single frame. This is a general problem with using numerical solutions. They are simply inaccurate. So what do we do?

In the case of the ball, we can avoid using a numerical solution altogether, and go to an analytical solution. I won't go through how I got the solution, but here is the different section:

float h0;
float t = 0;
float pd;

void setup() {
  size(400, 400);
  iterator = 0;
  yHist = new float[width];
  noFill();
  h0 = height - y;
  pd = 2*sqrt(h0/g);
}

void draw() {
  background(255);
  y = g*sq((t-pd/2)%pd - pd/2) + height - h0;
  t += 0.5;
  ellipse(x, y, radius*2, radius*2);
... etc.

This is all well and good for a bouncing ball, but a double pendulum is a much more complex system. There is no fully analytical solution to the double pendulum problem. So how do we minimize error in a numerical solution?

One strategy is to take smaller steps. The smaller the steps you take, the closer you are to the real solution. You can do this by reducing g (this might feel like cheating, but think for a minute about the units you're using. g=9.81 m/s^2. How does that translate to pixels and frames?). This will also make the pendulum move slower on the screen. If you want to increase accuracy without changing the viewing pace, you can take many small steps before rendering the frame. Consider changing lines 39-46 to

int substepCount = 1000;
for (int i = 0; i < substepCount; i++) {
  acc1 = (g*(sin(ang2)*cos(ang1-ang2)-mu*sin(ang1))-(l2*vel2*vel2+l1*vel1*vel1*cos(ang1-ang2))*sin(ang1-ang2))/(l1*(mu-cos(ang1-ang2)*cos(ang1-ang2)));
  acc2 = (mu*g*(sin(ang1)*cos(ang1-ang2)-sin(ang2))+(mu*l1*vel1*vel1+l2*vel2*vel2*cos(ang1-ang2))*sin(ang1-ang2))/(l2*(mu-cos(ang1-ang2)*cos(ang1-ang2)));
    
  vel1 += acc1/substepCount;
  vel2 += acc2/substepCount;
  
  ang1 += vel1/substepCount;
  ang2 += vel2/substepCount;
}

This changes your one big step to 1000 smaller steps, making it much more accurate. I tested that part out and it continued for over 20000 frames multiple times with no NaN errors. It might devolve into NaN at some point, but this allows it to last much longer.

EDIT:

I also highly recommend using % TWO_PI when incrementing the angles:

ang1 = (ang1 + vel1/substepCount) % TWO_PI;
ang2 = (ang2 + vel2/substepCount) % TWO_PI;

because it makes the angle measurements MUCH more accurate in the later times.

When you don't do this, if vel1 is positive for a long time, then ang1 gets bigger and bigger. Once ang1 is greater than 1, the computer needs a bit to indicate the ones place, at the expense of an extra digit on the end. Since numbers are stored using binary, this happens again when ang1 > 2, and again when ang1 > 4, and so on.

If you keep it between -PI and PI (which is what % does in this case), you only need a bit for the sign and a bit for the ones place, and all the remaining bits can be used to measure the angle to the highest possible precision. This is actually important: if vel1/substepCount < 1/32768, and ang1 doesn't have enough bits to measure out to the 1/32768 place, then ang1 will not register the change.

To see the effects of this difference, give ang1 and ang2 really high initial values:

g = 0.0981;
ang1 = 101.1*PI;
ang2 = 101.1*PI;

If you don't use % TWO_PI, it approximates low velocities to zero, resulting in a bunch of stopping and starting.

END EDIT

If you need it to go for a ridiculously long time, so long that it isn't feasible to increase substepCount sufficiently, there is another thing you can do. This all comes about because vel increases to an extreme degree. You can constrain vel1 and vel2 so that they don't get too big.

In this case, I would recommend limiting the velocities based on conservation of energy. There is a maximum amount of mechanical energy allowed in the system based on the initial conditions. You cannot have more mechanical energy than the initial potential energy. Potential energy can be calculated based on the angles:

U(ang1, ang2) = -g*((m1+m2)*l1*cos(ang1) + m2*l2*cos(ang2))

Therefore we can determine exactly how much kinetic energy is in the system at any moment: The initial values of ang1 and ang2 give us the total mechanical energy. The current values of ang1 and ang2 give us the current potential energy. Then we can simply take the difference in order to find the current kinetic energy.

The way that pendulum motion is typically described does not lend itself to computing kinetic energy. It is possible, but I'm not going to do it here. My recommendation for constraining the velocities of the two pendulums is as follows:

1. Calculate the kinetic energy of the two arms separately.
2. Take the ratio between them
3. Calculate the total kinetic energy currently in the two arms.
4. Distribute the kinetic energy in the same proportions as you calculated in step 2. e.g. If you calculate that there is twice as much kinetic energy in the further mass as there is in the closer mass, put 1/3 of the kinetic energy in the closer mass and 2/3 in the further one.

I hope this helps, let me know if you have any questions.