如何计算摆的运动？

Question

我正在尝试使用 Runge-Kutta 四阶方法对钟摆进行简单的模拟。我正在使用 p5.js。通常它会正确计算角度，但有时它只是开始随机旋转等。我不知道如何解决这个问题，我认为这似乎是我实现上述算法的问题？

我使用了庞涛的计算物理学导论中的代码作为灵感，看起来很相似

上述代码在 Fortran 中的第一部分

第二部分

let screenWidth = 1300;
let screenHight = 1970;
let angleChangeDifference;

let gSlider;
let lSlider;
let aFrequencySlider;
let dumpingSlider;
let startAngleSlider;
let timestepSlider;
let timeMaxSlider;
let dForceSlider;
let initForceSlider;

let startAngle = 0;
let currentAngle = 0;

let circleX = 0;
let circleY = 200;
let circleRWidth = 100;
let circleRHeight = 100;

let lineXStart = 0;
let lineYStart = 0;
let lineXEnd = circleX;
let lineYEnd = (circleY - circleRHeight / 2);

function setup() {
  createCanvas(screenWidth, screenHight);
  changeScreenDeafultStartingPoint(screenWidth / 2, 100);
  frameRate(60)
  createSliders();
  setInterval(showSliderValue, 100);
  startButton.mouseClicked(start);
  restartButton.mouseClicked(restart);
  chartButton.mouseClicked(enableChart);
  background(200);
}
let angleSign = '\u00B0';
let omegaSign = '\u03C9';
let chartOn = false;
let step = 0;

function draw() {
  startupConfiguration()
  showSliderValue()
  if (step >= 1) {
    rotatePendulum();
    if (step == 1) {}
    if (step == 2)
      startButton.remove();
    chartButton.position(20, 340);
    if (chartOn == true) {
      createChart(0, 1, 'czas (s)', 'kat (' + angleSign + ')', degreesArr);
      createChart(0, 320, 'czas (s)', omegaSign + ' (' + angleSign + '/s)', omegaArr);
      scale(2);
      createPhaseChart(290, 80, 'kat (' + angleSign + ')', omegaSign + ' (' + angleSign + '/s)', degreesArr);
      scale(0.5);
    }
  }
  line(lineXStart, lineYStart, lineXEnd, lineYEnd + (20 * lSlider.value()));
  fill(200, 76, 43)
  ellipse(circleX, circleY + (20 * lSlider.value()), circleRWidth, circleRHeight);
}

function createSliders() {

  gSlider = createSlider(0.05, 20, 9.81, 0.01);
  gSlider.position(1100, -90);

  lSlider = createSlider(0.5, 10, 5, 0.5);
  lSlider.position(1100, -50);

  initForceSlider = createSlider(-5, 5, 0, 0.5);
  initForceSlider.position(1100, 50);

  dForceSlider = createSlider(-2, 2, 0.9, 0.05)
  dForceSlider.position(1100, 90);

  aFrequencySlider = createSlider(-2, 2, 2 / 3, 1 / 3);
  aFrequencySlider.position(1100, 130);

  dumpingSlider = createSlider(0.00, 1.5, 0.5, 0.05);
  dumpingSlider.position(1100, 170);

  startAngleSlider = createSlider(-Math.PI / 2, Math.PI / 2, 0, Math.PI / 32);
  startAngleSlider.position(1100, 210);

  timestepSlider = createSlider(0, 1000, 100, 10);
  timestepSlider.position(1100, 250);

  timeMaxSlider = createSlider(10, 10000, 1000, 10);
  timeMaxSlider.position(1100, 290);


  startButton = createButton('ZATWIERDZ', false);
  startButton.position(100, 310)

  restartButton = createButton('RESTART', false);
  restartButton.position(20, 310)

  chartButton = createButton('WYKRES', false);
  chartButton.position(-200, 340);

}

function showSliderValue() {
  background(200);
  fill(0, 0, 0)

  text('sila poczatkowa', 440, -60)
  text(initForceSlider.value(), 400, -42)
  text('sila sprawcza', 440, -20)
  text(dForceSlider.value(), 400, -2)
  text('czestosc katowa', 440, 20)
  text(aFrequencySlider.value(), 400, 42)

  text('tlumienie', 440, 60)
  text(dumpingSlider.value(), 400, 82)
  text('kat poczatkowy', 440, 100)
  text(int(degrees(startAngleSlider.value())), 400, 122)
  text('krok czasowy (N1)', 440, 140)
  text(timestepSlider.value(), 400, 162)
  text('dlugosc symulacji (N2)', 440, 180)
  text(timeMaxSlider.value(), 400, 202)


}

function start() {
  angleIndex = 0;
  step++;
  startAngle = startAngleSlider.value();
  currentAngle = startAngle;
  angleChangeDifference = simulate();
  rotatePendulum()

  startButton.html("START")
}

function restart() {
  window.location.reload();
}



function enableChart() {
  chartOn = true;
}

function createChart(moveByX, moveByY, xName, yName, table) {
  rotate(-currentAngle);
  scale(1.1);
  translate(moveByX, moveByY);
  strokeWeight(1);
  line(-500, 500, 530, 500);
  line(-500, 700, 530, 700);
  line(-500, 500, -500, 700);
  line(530, 500, 530, 700);
  strokeWeight(1);
  let counter = 0;
  for (i = 510; i < 700; i += 10) {
    if (counter < 9 && counter % 2 == 0)
      text(90 - 10 * counter, -520, i + 5)
    else if (counter == 9 && counter % 2 == 0)
      text(90 - 10 * counter, -515, i + 5)
    else if (counter > 9 && counter % 2 == 0)
      text(90 - 10 * counter, -525, i + 5)

    line(-505, i, 530, i);
    counter++;
  }
  textSize(25);
  text(xName, -20, 750)
  textSize(12);

  counter = 0;
  for (i = -490; i < 535; i += 25) {
    line(i, 500, i, 705);
    if (counter % 4 == 0) {
      line(i, 500, i, 705);
      text(counter * 2.5, i - 5, 715);
    }
    counter++;
  }
  rotate(-90);
  textSize(25);
  text(yName, -670, -550)
  textSize(12);
  rotate(90);

  fillChartByTableValues(table);
  translate(-moveByX, -moveByY);
  scale(0.91);
  rotate(currentAngle);
}

function fillChartByTableValues(table) {
  strokeWeight(2);
  stroke(0, 0, 255);
  for (i = 0; i < timeArr.length - 1; i++) {
    FirstPointX = -490 + timeArr[i] * 10;
    FirstPointY = 600 + table[i] * (-1);
    SecondPointX = -490 + timeArr[i + 1] * 10;
    SecondPointY = 600 + table[i + 1] * (-1);
    line(FirstPointX, FirstPointY, SecondPointX, SecondPointY);
  }
  stroke(0, 0, 0);
  strokeWeight(0.1);
}


function createPhaseChart(moveByX, moveByY, xName, yName, table) {
  rotate(-currentAngle);
  scale(1.1);
  translate(moveByX, moveByY);
  strokeWeight(1);
  line(-500, 500, -300, 500);
  line(-500, 700, -300, 700);
  line(-500, 500, -500, 700);
  line(-300, 500, -300, 700);
  strokeWeight(1);
  let counter = 0;
  textSize(8);
  for (i = 510; i < 700; i += 10) {
    if (counter < 9 && counter % 2 == 0)
      text(90 - 10 * counter, -517, i + 3)
    else if (counter == 9 && counter % 2 == 0)
      text(90 - 10 * counter, -512, i + 3)
    else if (counter > 9 && counter % 2 == 0)
      text(90 - 10 * counter, -520, i + 3)

    line(-505, i, -300, i);
    counter++;
  }
  textSize(12);

  textSize(15);
  text(xName, -430, 735)
  textSize(12);

  counter = 0;
  textSize(8);
  for (i = -490; i < -300; i += 10) {
    line(i, 500, i, 705);

    if (counter < 9 && counter % 2 == 0)
      text(-90 + 10 * counter, i - 7, 715);
    else if (counter == 9 && counter % 2 == 0)
      text(-90 + 10 * counter, i - 2, 715);
    else if (counter > 9 && counter % 2 == 0)
      text(-90 + 10 * counter, i - 4, 715);

    counter++;
  }
  textSize(12);

  rotate(-90);
  textSize(15);
  text(yName, -620, -528)
  textSize(12);
  rotate(90);

  fillPhaseChartByTableValues(degreesArr, omegaArr)

  translate(-moveByX, -moveByY);
  scale(0.91);
  rotate(currentAngle);
}


function fillPhaseChartByTableValues(tableX, tableY) {
  translate(-400, 600);
  strokeWeight(1);
  stroke(0, 0, 255);
  for (i = 0; i < tableX.length; i++) {
    ellipse(tableX[i], tableY[i], 0.5, 0.5);
  }
  translate(400, -600);
  stroke(0, 0, 0);
  strokeWeight(0.1);
}


function startupConfiguration() {
  background(200);
  angleMode(DEGREES);
  changeScreenDeafultStartingPoint(screenWidth / 2, 100);
}

function changeScreenDeafultStartingPoint(x, y) {
  translate(x, y);
}


let angleIndex = 0;

function rotatePendulum() {
  currentAngle = angleChangeDifference[angleIndex] * (180 / PI);
  rotate(currentAngle);
  if (step > 1) {
    angleIndex++
  }
}

function calculateIntegral(t, q, dt, f) {

  let k1 = f(t, q).map(val => val * dt);

  let temp = k1.map(val => val * 0.5);
  temp = temp.map((val, index) => val + q[index])

  let k2 = f(t + 0.5 * dt, temp).map(val => val * dt);

  temp = k2.map(val => val * 0.5);
  temp = temp.map((val, index) => val + q[index])
  let k3 = f(t + 0.5 * dt, temp).map(val => val * dt);

  temp = q.map((val, index) => val + k3[index]);
  let k4 = f(t + dt, temp).map(val => val * dt);

  temp = k2.map((val, index) => val + k3[index])
  temp = temp.map(val => val * 2)
  temp = temp.map((val, index) => (val + k1[index] + k4[index]) / 6)

  temp = temp.map((val, index) => val + q[index])
  return [t + dt, temp];
}


function modelPendulum(t, q) {
  let c = dumpingSlider.value();
  let fw = dForceSlider.value();
  let w = aFrequencySlider.value();

  let x1 = q[0];
  let x2 = q[1];
  return [x2, -(Math.sin(x1)) - (c * x2) + (fw * Math.cos(w * t))];
}

let degreesArr, timeArr;

function simulate() {
  let t = 0.0;
  let dt = (3 * Math.PI) / timestepSlider.value()
  let tf = timeMaxSlider.value() * dt
  let q = [startAngle, initForceSlider.value()];
  let Nt = int(Math.round((tf - t) / dt)) + 1;
  let solution = new Array(q.length + 1);

  for (i = 0; i < q.length + 1; i++) {
    solution[i] = new Array(Nt).fill(0);
  }
  solution[0][0] = t;
  solution[1][0] = q[0];
  solution[2][0] = q[1];

  k = 1;
  while (t <= tf) {
    let temporaryResult = [];
    temporaryResult = calculateIntegral(t, q, dt, modelPendulum);
    t = temporaryResult[0];
    q = temporaryResult[1];
    solution[0][k] = t;
    solution[1][k] = q[0];
    solution[2][k] = q[1];

    k = k + 1
  }

  timeArr = solution[0];
  degreesArr = solution[1];
  omegaArr = solution[2];
  let counter = 0;
  let ifChaos = false;
  while (counter != degreesArr.length - 1 && ifChaos != true) {
    if (degreesArr[counter] > 13.5 || degreesArr[counter] < -13.5) {
      ifChaos = true;
    }
    counter++;
  }

  if (ifChaos == true) {
    degreesArr = degreesArr.map(val => val * 5.32);
    omegaArr = omegaArr.map(val => val * 23.32);
  } else {
    degreesArr = degreesArr.map(val => val * 35.32);
    omegaArr = omegaArr.map(val => val * 35.32);
  }
  return solution[1]
}

<script src="https://cdn.jsdelivr.net/npm/p5@1.1.9/lib/p5.js"></script>

Answer 1

我对 P5 不熟悉，所以这里有一个使用 ThreeJS 的解决方案，利用来自 Mathematics Stack Exchange 的 Runge-Kutta 算法。

为方便起见，我将 Runge-Kutta 算法封装在一个类中，构造函数采用以下初始参数：

• 重力加速度常数g（对于地球，9.81 米/秒/秒），
• 摆锤长度（以米为单位），
• 摆锤的初始角度（0 为直下），
• 初始角速度（以米/秒为单位），以及
• 最大时间增量。 （由于 Runge-Kutta 用于求解以时间为变量的二阶微分方程，因此根据作者的实验，似乎不能过度扩展 delta 时间增量并仍然保持所得钟摆位置和速度的准确性。此参数只是限制updatePosition方法中传递的最大t值，默认为0.1s。）

为了帮助使用 Runta-Kutta 算法，下面的代码模拟了两个 1 米钟摆：

• 第一个初始位置为 -PI/2（-90 度），没有起始角速度。
• 第二个初始位置为 PI (180 度)，起始角速度非常小。

<script type="module">

  import * as THREE from 'https://cdn.jsdelivr.net/npm/three@0.115.0/build/three.module.js';

  class RungeKutta {
  
    constructor( g, pendulumLength, initialAngle, angularVelocity, maxTimeDelta ) {
    
      this.g = g;
      this.pendulumLength = pendulumLength;
      this.theta = initialAngle;
      this.omega = angularVelocity;
      this.maxTimeDelta = maxTimeDelta || 0.1;
      
    }
    
    updatePosition( t ) {
    
      let self = this;
    
      function omegaDot( theta ){
        return -( self.g / self.pendulumLength ) * Math.sin( theta );
      }

      function thetaDot( omega ){
        return omega;
      }   

      // If the browser tab becomes inactive, then there will be a large
      // time delta, which will disrupt the RungeKutta algorithm.  If more
      // than max allowed seconds has lapsed, then reset the timer.
      if ( self.maxTimeDelta < t ) {
        t = self.maxTimeDelta;
      }
      
      let aomega = omegaDot( self.theta );
      let atheta = thetaDot( self.omega );
      let bomega = omegaDot( self.theta + 0.5 * t * atheta );
      let btheta = thetaDot( self.omega + 0.5 * t * aomega );
      let comega = omegaDot( self.theta + 0.5 * t * btheta );
      let ctheta = thetaDot( self.omega + 0.5 * t * bomega );
      let domega = omegaDot( self.theta + t * ctheta );
      let dtheta = thetaDot( self.omega + t * comega );

      self.omega = self.omega + ( t / 6 ) * ( aomega + 2 * bomega + 2 * comega + domega );
      self.theta = self.theta + ( t / 6 ) * ( atheta + 2 * btheta + 2 * ctheta + dtheta );

      return self;
      
    }
    
  }

  //
  // Set up the ThreeJS environment.
  //
  var renderer = new THREE.WebGLRenderer();
  renderer.setSize( window.innerWidth, window.innerHeight );
  document.body.appendChild( renderer.domElement );

  var camera = new THREE.PerspectiveCamera( 45, window.innerWidth / window.innerHeight, 1, 500 );
  camera.position.set( 0, 0, 100 );
  camera.lookAt( 0, 0, 0 );

  var scene = new THREE.Scene();

  //
  // Create the pendulum mesh.
  //
  var length = 30, width = 1;

  var shape = new THREE.Shape();
  shape.moveTo( -width / 2, 0 );
  shape.lineTo( +width / 2, 0 );
  shape.lineTo( +width / 2, -length );
  shape.lineTo( -width / 2, -length );
  shape.lineTo( -width / 2, 0 );

  var extrudeSettings = {
    steps: 2,
    depth: 2,
    bevelEnabled: true,
    bevelThickness: .25,
    bevelSize: .25,
    bevelOffset: 0,
    bevelSegments: 1
  };

  var geometry = new THREE.ExtrudeBufferGeometry( shape, extrudeSettings );
  var material = new THREE.MeshBasicMaterial( { color: 0x00ff00 } );
  var mesh0 = new THREE.Mesh( geometry, material );
  mesh0.position.x = -17;
  scene.add( mesh0 );
  var mesh1 = mesh0.clone();
  mesh1.position.x = +17;
  scene.add( mesh1 );
  
  //
  // And now animate the pendulum using RungeKutta.
  //
  let pendulumState0 = new RungeKutta( 9.81, 1, -Math.PI / 2, 0,    0.1 );
  let pendulumState1 = new RungeKutta( 9.81, 1, Math.PI,      0.01, 0.1 );
  
  let now = performance.now();
  let lastTimer = now;
  
  var animate = function () {
  
    requestAnimationFrame( animate );
    
    now = performance.now();
    pendulumState0.updatePosition( ( now - lastTimer ) / 1000 );
    pendulumState1.updatePosition( ( now - lastTimer ) / 1000 );
    lastTimer = now;

    mesh0.rotation.z = pendulumState0.theta;
    mesh1.rotation.z = pendulumState1.theta;
    
    renderer.render( scene, camera );
    
  };

  animate();
</script>

希望这将有助于您的 P5 实施。

Answer 2

这是一个带有 p5 的 pendulum 的示例，使用 Leapfrog method，但它可以很容易地适应使用 Runge-Kutta 4（Leapfrog 的优势在于它是辛的，与 RK 不同）。

对于一个简单的钟摆，您必须求解二阶微分（牛顿或欧拉-拉格朗日）方程

d^2 O / dt^2 = -(g/l) sin(O)。

其中 O 是与垂直方向的夹角。

应用数值方法首先方便的是把它转化为两个一阶微分方程的等价系统

dw/dt = -(g/l) sin(O)

dO/dt = w

您可以使用 Runge-Kutta4 或其他方法解决此系统。一个非常简单的方法是简单的二阶 Leapfrog 方法。你在步骤 dt 中离散时间所以

w_n == w((n-1/2)dt)

O_n == O(n dt)

然后在循环中迭代地执行以下操作

w_{n+1} = w_n -(g/l) sin(O_n) dt

O_{n+1} = O_n + w_{n+1} dt

这个方法在这个code的一个摆类中实现，连同它的参数，一个非常简单的render()函数，用于在p5中进行可视化，以及一个计算能量的方法，例如（可用于检查守恒能量，如果不强迫）。它非常简单，您可以对其进行测试。当然，应该使用小的时间步长 dt 以获得良好的准确性。

如果您希望使用 Runge-Kutta4 或其他高阶方法，可以方便地编写函数 Fw(t,w,O) 和 FO(t,w,O) 以便读取方程

dw/dt = Fw(t,w,O) (== -(g/l) sin(O) + F0 cos(W t))

dO/dt = FO(t,w,O) (== w )

我现在在其中包含了一个强制 F0 cos(W t)，以显示如何添加任何时间相关的力。

然后您可以及时离散化（没有越级步骤）

t_n == n dt

w_n == w(n dt)

O_n == O(n dt)

并计算每个时间步 n 的RK4 的数量 k1、k2、k3、k4。 请注意，对于这种情况，ki 是二维向量 {kiw,kiO}，函数 f(t,w,O) =={Fw(t,w,O),FO(t,w,O)} 是也是一个二维向量。

考虑到这一点，并分别定义函数 Fw 和 FO，您可以轻松地将 p5 示例 1 的 Leapfrog 方法替换为 RK4 方法，并在需要时添加强制。

如果您看到非常奇怪的行为，您应该检查标志并正确设置 RK4 方法。

我希望这会有所帮助。

let P1;

function setup() {
  createCanvas(720, 400);
  P1=new PenduloSimple(3.14,0.0);
}

function draw() {
  background(220);
  P1.render();
  P1.leapFrog();
}


class PenduloSimple {
  constructor(ang, velang) {
    this.ang = ang;
    this.velang = velang;
    this.g = 9.8;
    this.dt = 0.01;
    this.l = 1.0;
    this.m = 0.1;
    this.E0 = this.Energy();
  }
  Energy() {
    let E=
    this.m * this.l * this.l * this.velang * this.velang * 0.5 - this.m * this.g * this.l * cos(this.ang);
    return E;    
  }
  leapFrog() {
    // Método de Leapfrog
    this.velang = 
    this.velang + this.dt*(-this.g/this.l)*sin(this.ang);
    this.ang = this.ang + this.dt*this.velang;
  }  
  render(){
    var x0 = width / 2;
    var y0 = height / 2;
    var mult0 = 100;  

    var x = this.l * sin(this.ang);
    var y = this.l * cos(this.ang);

    var xplot = x0 + x * mult0;
    var yplot = y0 + y * mult0;

    // Draw a circle
    stroke(50);
    fill(100);
    ellipse(xplot, yplot, 24, 24);

    // La cuerda
    stroke(50);
    line(x0, y0, xplot, yplot);

    let E = this.Energy();

    stroke(0);
    fill(255, 0, 0);
    rect(10, height * 0.5, 20, -E * height * 0.45 / this.E0);
    fill(0, 255, 0);
    rect(30, height * 0.5, 20, -this.ang * height * 0.4 / 6.28);
    fill(0, 0, 255);
    rect(50, height * 0.5, 20, -this.velang * mult0 * 0.1);

    fill(1);
    stroke(255);
    fill(255, 0, 0);
    text("E: " + E, width * 0.1, 20);
    fill(50, 200, 50);
    text("ang: " + this.ang, width * 0.1, 30);
    fill(0, 0, 255);
    text("velang: " + this.velang, width * 0.1, 40);
  }
}