Force calculation for magnetic pendulum and numerical integration

Hello everyone,
I’ve been working on a simulation of a magnetic pendulum (Magnetic Pendulum for reference). Now, I have a question concerning how you calculate the forces/acceleration:
In examples you find online of the magnetic pendulum (such as the one provided above), and other physics force and acceleration calculations, the forces are calculated at a starting position, the position updated according to time-step dt and then new forces at time t+dt calculated to get a new position and so on. Basic numerical integration, makes sense.

However, I noticed, that if I calculate the forces this way, things don’t go as one would expect. For example: when the pendulum is attracted by a magnet it just stops at the particluar magnet even tough you would expect it to “overswing” a bit. Therefore, I introduced three variables for each force which is all previous forces added together, and now it somehow seems to work correctly.
I’m wondering if

  1. the way of calculating the forces I used makes sense in a physics simulation and
  2. if not, then what am I doing wrong? Because it seems to work for other people.

Here’s my code:

sketch.js
p5.disableFriendlyErrors = true;
let magnete = [];
let particles = [];
var maganzahl = 3; //number of magnets
var magradius = 100; //radius of magnets from mid-point
var G = 0.002; //coefficient of force towards middle (gravity)
var R = 0.2; //friction coefficient
var h = 50; //height of bob over magnets
var M = 150000; //strength of magnets
var m = 1; //mass (might not work)
var dt = 0.25; //timestep
var d = 2; //pixel density
var counter2 = 0;
var Reset;
var end = false;
//--------------------------------------------

function setup() {
  createCanvas(600, 600);
  pixelDensity(d);
  background(60);

  Reset = createButton('Reset');
  Reset.position(width + 19, 49);
  Reset.mousePressed(resetcanvas);

  //construction of magnets
  for (var i = 0; i < maganzahl; i++) {
    magnete.push(new Magnet((width / 2) + magradius * cos(TWO_PI * (i / maganzahl)), (height / 2) + magradius * sin(TWO_PI * (i / maganzahl)), i));
  }
}

//construction of a new "starting position" by mouse-click
function mousePressed() {
  if (mouseX < width && mouseY < height) {
    particles.push(new Particle(mouseX, mouseY));
  }
}

function draw() {

  for (var i = 0; i < particles.length; i++) {
    if (particles[i].counter < 1000) {
      //5 updates per frame(to speed it up)
      for (var k = 0; k < 5; k++) {
        for (var j = 0; j < magnete.length; j++) {
          particles[i].attracted(magnete[j]);
        }
        particles[i].update();
        particles[i].show2();
      }
    } else if (particles[i].counter < 1001) {
      particles[i].counter += 1;

      var nearest = [];
      for (var j = 0; j < magnete.length; j++) {
        nearest.push(particles[i].near(magnete[j]));
      }
      if (nearest.indexOf(min(nearest)) == 0) {
        var c = color("green");
      }
      if (nearest.indexOf(min(nearest)) == 1) {
        var c = color("purple");
      }
      if (nearest.indexOf(min(nearest)) == 2) {
        var c = color("orange");
      }
      if (nearest.indexOf(min(nearest)) == 3) {
        var c = color("blue");
      }
      if (nearest.indexOf(min(nearest)) == 4) {
        var c = color("red");
      }
      if (nearest.indexOf(min(nearest)) == 5) {
        var c = color("yellow");
      }
      //show particle trace according to nearest magnet
      particles[i].show(c);
    }
  }
  //displaying magnets
  for (var i = 0; i < magnete.length; i++) {
    magnete[i].show();
  }
  //displaying mid-point
  stroke(255);
  circle(width / 2, height / 2, 3);
}

function resetcanvas() {
  background(60);
}p5.disableFriendlyErrors = true;
let magnete = [];
let particles = [];
var maganzahl = 3; //number of magnets
var magradius = 100; //radius of magnets from mid-point
var G = 0.002; //coefficient of force towards middle (gravity)
var R = 0.2; //friction coefficient
var h = 50; //height of bob over magnets
var M = 150000; //strength of magnets
var m = 1; //mass (might not work)
var dt = 0.25; //timestep
var d = 2; //pixel density
var counter2 = 0;
var Reset;
var end = false;
//--------------------------------------------

function setup() {
  createCanvas(600, 600);
  pixelDensity(d);
  background(60);

  Reset = createButton('Reset');
  Reset.position(width + 19, 49);
  Reset.mousePressed(resetcanvas);

  //construction of magnets
  for (var i = 0; i < maganzahl; i++) {
    magnete.push(new Magnet((width / 2) + magradius * cos(TWO_PI * (i / maganzahl)), (height / 2) + magradius * sin(TWO_PI * (i / maganzahl)), i));
  }
}

//construction of a new "starting position" by mouse-click
function mousePressed() {
  if (mouseX < width && mouseY < height) {
    particles.push(new Particle(mouseX, mouseY));
  }
}

function draw() {

  for (var i = 0; i < particles.length; i++) {
    if (particles[i].counter < 1000) {
      //5 updates per frame(to speed it up)
      for (var k = 0; k < 5; k++) {
        for (var j = 0; j < magnete.length; j++) {
          particles[i].attracted(magnete[j]);
        }
        particles[i].update();
        particles[i].show2();
      }
    } else if (particles[i].counter < 1001) {
      particles[i].counter += 1;

      var nearest = [];
      for (var j = 0; j < magnete.length; j++) {
        nearest.push(particles[i].near(magnete[j]));
      }
      if (nearest.indexOf(min(nearest)) == 0) {
        var c = color("green");
      }
      if (nearest.indexOf(min(nearest)) == 1) {
        var c = color("purple");
      }
      if (nearest.indexOf(min(nearest)) == 2) {
        var c = color("orange");
      }
      if (nearest.indexOf(min(nearest)) == 3) {
        var c = color("blue");
      }
      if (nearest.indexOf(min(nearest)) == 4) {
        var c = color("red");
      }
      if (nearest.indexOf(min(nearest)) == 5) {
        var c = color("yellow");
      }
      //show particle trace according to nearest magnet
      particles[i].show(c);
    }
  }
  //displaying magnets
  for (var i = 0; i < magnete.length; i++) {
    magnete[i].show();
  }
  //displaying mid-point
  stroke(255);
  circle(width / 2, height / 2, 3);
}

function resetcanvas() {
  background(60);
}
particles.js
function Particle(x, y) {
  this.orgpos = createVector(x, y);
  this.pos = createVector(x, y);
  this.prev = createVector(x, y);
  this.vel = createVector();
  this.acc = createVector();
  this.accpre = createVector();
  this.accprepre = createVector();
  this.velprediction = this.vel;
  this.magnetf = createVector();
  this.gravity = createVector();
  this.friction = createVector();
  this.shape = new Array();
  this.counter = 0;


  //calculating new positions
  this.update = function() {
    //predictor for velocity -> Beeman's algorithm
    this.velprediction.add(this.accpre.mult(3 / 2 * dt).add(this.accprepre.mult(-1 / 2 * dt)));

    var momgrav = createVector(width / 2 - this.pos.x, height / 2 - this.pos.y);
    var momfric = createVector(this.velprediction.x, this.velprediction.y);
    momgrav.mult(G * m); //force due to gravity
    momfric.mult(-R); //force due to friction
    this.gravity.add(momgrav);
    this.friction.add(momfric);

    //a = F/m
    this.acc = createVector((this.magnetf.x + this.gravity.x + this.friction.x) / m, (this.magnetf.y + this.gravity.y + this.friction.y) / m);

    //-=Beeman's Algorithm=-
    this.vel.add(this.acc.mult(dt * 1 / 3).add(this.accpre.mult(dt * 5 / 6)).add(this.accprepre.mult(-1 / 6 * dt)));
    this.pos.add(this.vel.mult(dt).add(this.accpre.mult(dt * dt * 2 / 3)).add(this.accprepre.mult(-1 / 6 * dt * dt)));

    this.accprepre = createVector(this.accpre.x, this.accpre.y);
    this.accpre = createVector(this.acc.x, this.acc.y);
    this.velprediction = createVector(this.vel.x, this.vel.y);
    this.counter += 1;
    this.shape.push(new p5.Vector(this.pos.x, this.pos.y));
  }

  //calculating force due to magnets -> attracted called earlier than update in sketch.js
  this.attracted = function(target) {
    var magn = createVector(target.pos.x - this.pos.x, target.pos.y - this.pos.y);
    var dist = sqrt(sq(h) + sq(magn.x) + sq(magn.y)); //distance bob - magnet
    strength = M / (Math.pow(dist, 3));
    magn.mult(strength);
    this.magnetf.add(magn);
  }

  //calculating distance to target
  this.near = function(target) {
    var dist = sqrt(sq(h) + sq(this.pos.x - target.pos.x) + sq(this.pos.y - target.pos.y));
    return (dist);
  }

  //trace
  this.show = function(col) {
    beginShape();
    stroke(col);
    for (var i = 0; i < this.shape.length - 1; i++) {
      line(this.shape[i].x, this.shape[i].y, this.shape[i + 1].x, this.shape[i + 1].y);
      strokeWeight(2);
      noFill();
    }
    endShape();
  }
  //dots
  this.show2 = function() {
    strokeWeight(1)
    point(this.pos.x, this.pos.y);
  }
}
magnete.js
function Magnet(x, y, n) {
  this.pos = createVector(x, y);
  this.n = n;

  this.show = function() {
    noStroke();
    //color for each magnet
    if (n == 0) {
      fill("green");
    }
    if (n == 1) {
      fill("purple");
    }
    if (n == 2) {
      fill("orange");
    }
    if (n == 3) {
      fill("blue");
    }
    if (n == 4) {
      fill("red");
    }
    if (n == 5) {
      fill("yellow");
    }
    strokeWeight(4);
    circle(this.pos.x, this.pos.y, 10);

  }
}

p5.js Web Editor

Any help would be greatly appreciated!

1 Like