为什么我的双摆代码返回 NaN？

如何解决为什么我的双摆代码返回 NaN？

当我打印我的加速度和速度时，它们都开始（看似）正常。不久，它们开始变得非常大，然后返回 -Infinity，然后返回 NaN。我在数学/物理方面已经尽力了，但我的知识有限，所以要温柔。任何帮助将不胜感激。

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,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); // circle at m2
  
  line(0,l1*cos(ang1),l1*sin(ang1)); // arm 1
  line(l1*cos(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,vel2);
  
  popMatrix();
}

解决方法

您的代码没有做错任何事情，但是这种数学技术的应用很棘手。

这是对微分方程使用数值“解”的一般问题。如果您尝试模拟弹跳球，则会发生类似的事情：

//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;
}

如果您运行该代码，您会注意到球似乎每次都弹得更高。显然，这不会发生在现实生活中，即使在理想的无损场景中也不应该发生。那么这里发生了什么？

如果您曾经使用过 Euler's method 来求解微分方程，那么您可能会对这里发生的事情有所了解。实际上，当我们编写微分方程的模拟代码时，我们正在应用欧拉方法。在弹跳球的情况下，实际曲线向下凹（除了弹跳时的点）。当实解向下凹时，欧拉方法总是给出高估计。这意味着每一帧，计算机都会猜测一点太高。这些错误加起来，球弹得越来越高。

同样，对于你的钟摆，几乎每一帧都会获得一点更多的能量。这是使用数值解的普遍问题。他们只是不准确。那我们该怎么办？

在球的情况下，我们可以完全避免使用数值解，而转而使用解析解。我不会详细介绍我是如何得到解决方案的，但这里是不同的部分：

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,radius*2);
... etc.

对于弹跳球来说这一切都很好，但双摆是一个复杂得多的系统。双摆问题没有完全解析解。那么我们如何最小化数值解中的误差？

一种策略是采取更小的步骤。您采取的步骤越小，您就越接近真正的解决方案。您可以通过减少 g 来做到这一点（这可能感觉像是在作弊，但请考虑一下您正在使用的单位。g=9.81 m/s^2。这如何转换为像素和帧？）。这也会使钟摆在屏幕上的移动速度变慢。如果您想在不改变观看速度的情况下提高准确性，您可以在渲染帧之前采取许多小步骤。考虑将第 39-46 行更改为

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;
}

这会将您的一大步更改为 1000 个较小的步，从而使其更加准确。我测试了那部分，它持续了超过 20000 帧多次，没有出现 NaN 错误。它可能会在某个时候转变为 NaN，但这可以让它持续更长时间。

编辑：

我也强烈推荐在增加角度时使用 % TWO_PI：

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

因为它使以后的角度测量更加准确。

当你不这样做时，如果vel1长期为正值，那么ang1会越来越大。一旦 ang1 大于 1，计算机需要一点来指示个位，代价是末尾多出一个数字。由于数字使用二进制存储，因此在 ang1 > 2 时会再次发生这种情况，在 ang1 > 4 时会再次发生这种情况，依此类推。

如果你把它放在 -PI 和 PI 之间（这就是 % 在这种情况下所做的），你只需要一点点作为符号，一点点作为个位，并且所有剩余的位都可用于以尽可能高的精度测量角度。这实际上很重要：如果 vel1/substepCount < 1/32768 和 ang1 没有足够的位来测量到 1/32768 位，那么 ang1 将不会记录更改。

要查看这种差异的影响，请给 ang1 和 ang2 提供非常高的初始值：

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

如果不使用 % TWO_PI，它会将低速逼近为零，从而导致一堆停止和启动。

结束编辑

如果您需要它持续可笑很长时间，以至于无法充分增加 substepCount，您还可以做另一件事。这一切都是因为 vel 增加到一个极端程度。您可以限制 vel1 和 vel2，以免它们变得太大。

在这种情况下，我建议根据能量守恒限制速度。根据初始条件，系统中允许存在最大机械能。您不能拥有比初始势能更多的机械能。可以根据角度计算势能：

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

因此我们可以在任何时候准确地确定系统中有多少动能：ang1 和 ang2 的初始值为我们提供了总机械能。 ang1 和 ang2 的当前值为我们提供了当前的势能。然后我们可以简单地取差值来找到当前的动能。

钟摆运动的典型描述方式不适用于计算动能。这是可能的，但我不会在这里做。我对约束两个摆的速度的建议如下：

1. 分别计算两个臂的动能。
2. 取它们之间的比例
3. 计算当前两臂的总动能。
4. 按照您在步骤 2 中计算的相同比例分配动能。例如如果您计算出较远质量中的动能是较近质量中的动能的两倍，请将 1/3 的动能放在较近的质量中，将 2/3 放在较远的质量中。

希望对您有所帮助，如果您有任何问题，请告诉我。