How to calculate movement of pendulum?

Question

I am trying to make a simple simulation of pendulum using Runge–Kutta fourth-order method. I am using p5.js. Usually it calculates the angles properly but sometimes it just starts spinning randomly etc. I have no idea how to resolve this issue it seems to be a problem with my implementation of said algorithm I think?

I used code from pang tao's introduction to Computational Physics as an inspiration and it seems quite similar

First part of said code in Fortran

Second part

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/[email protected]/lib/p5.js"></script>

Answer 1

Am not familiar with P5, so here's a solution using ThreeJS leveraging the Runge-Kutta algorithm from the Mathematics Stack Exchange.

For convenience, I have wrapped the Runge-Kutta algorithm in a class, with the constructor taking the initial parameters of:

• the gravity acceleration constant g (for earth, 9.81 meters/sec/sec),
• the pendulum length (in meters),
• the initial angle of the pendulum (where 0 is straight down),
• the initial angular velocity (in meters/sec), and
• the max time increment. (Since the Runge-Kutta is employed to solve a second order differential equation using time as the variable, it appears based on this author's experimentation that one cannot overextend the delta time increment and still retain accuracy of the resulting pendulum position and velocity. This parameter simply limits the maximum t value passed in the updatePosition method, with a default of 0.1s.)

To assist in the use of the Runta-Kutta algorithm, the code below simulates two 1 meter pendulums:

• The first having an initial position of -PI/2 (-90 deg) with no starting angular velocity.
• The second having an initial position of PI (180 deg) with a very small starting angular velocity.

<script type="module">

  import * as THREE from 'https://cdn.jsdelivr.net/npm/[email protected]/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>

Hopefully this will assist with your P5 implementation.

Answer 2

This is an example of a pendulum with p5, using Leapfrog method, but it can be easily adapted to use Runge-Kutta 4 (Leapfrog has the advantage that it is symplectic, at variance with RK).

For a simple pendulum you have to solve the second order differential (Newton or Euler-Lagrange) equation

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

where O is the angle from the vertical.

To apply numerical methods is convenient first to translate it into an equivalent system of two first order differential equations

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

dO/dt = w

you can solve this system for instance with Runge-Kutta4 or other methods. A very simple one to start with is the simple 2nd order Leapfrog method. You discretize times in steps dt so

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

O_n == O(n dt)

and then do the following iteratively in a loop

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

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

This method is implemented in a pendulum class in this code, together with its parameters, a very simple render() function for visualization in p5, and a method for calculating energy for instance (can be used to check conservation of energy if not forced). It is very simple and you can test it. Of course, small time-steps dt should be used in order to have good accuracy.

If you wish to use Runge-Kutta4 or other high order method it is convenient to write functions Fw(t,w,O) and FO(t,w,O) such that the equations read

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

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

where I have now included a forcing F0 cos(W t) to show, just in case, how to add any time dependent force.

You can then discretize in time (no leapfrog steps)

t_n == n dt

w_n == w(n dt)

O_n == O(n dt)

and calculate the quantities k1, k2, k3, k4 of RK4 each time step n. Note that for this case the ki are two dimensional vectors {kiw,kiO}, and the function f(t,w,O) =={Fw(t,w,O),FO(t,w,O)} is also a 2d vector.

Taking this into account, and defining separately the functions Fw and FO you can easily replace the Leapfrog method of the p5 example 1 by a RK4 method, and add forcing if desired.

If you see very strange behaviour, you should check signs and a proper setup of the RK4 method.

I hope this can help.

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