Okay

In an effort to ask a more pointed question...

I've updated the VerletMover constructor to include a revised calculation for the initial state of the previousLocation variable. I've also made the code below more consistent and changed the display to graph the displacement value of each mover as a different coloured point over time.

My assumption was that the initial condition for the previousLocation variable could be appropriately set to (0,0) as no acceleration had been applied in the previous frame to the starting one. But even allowing for a previous acceleration calculation, the results are less appropriate not more so.

If you adjust the deltaTime variable as suggested in the comments you can see that the VerletMover class is the least stable over time while the EulerMover class tracks very closely to the NewtonianMover class. This is completely counter to my expectations.

So, to my questions(s). Have I misunderstood something in the implementation of the VerletMover class, perhaps the previousLocation state, that would account for this discrepancy? And is there a reason why the Euler class isn't deviating more?

Any help greatlly appreciated

Thanks

Matt

1. int margin;
   float deltaTime;
   float cumulativeTime;
   PVector globalAcceleration;
   NewtonianMover newtonianMover;
   EulerMover eulerMover;
   VerletMover verletMover;
   
   void setup() {
     size(1000, 600);
     smooth(5);
     background(0);
     frameRate(200);
   
     margin = 50;
   
     deltaTime = 1;//minimal deviation at this timestep but try 0.1 or 0.01
     cumulativeTime = 0;
     globalAcceleration = new PVector(10, 0);
   
     newtonianMover = new NewtonianMover();
     eulerMover = new EulerMover();
     verletMover = new VerletMover();
   }
   
   void draw() {
     translate(margin, height-margin);
     blendMode(ADD);
     cumulativeTime += deltaTime;
   
     newtonianMover.update();
     eulerMover.update();
     verletMover.update();
   
     newtonianMover.display();
     eulerMover.display();
     verletMover.display();
   
     if (frameCount%100==0) {
       println("Newtonian = "+newtonianMover.location.x+", Euler = "+eulerMover.location.x+", Verlet = "+verletMover.location.x);
     }
   }
   
   class NewtonianMover {//green line
     PVector location;
     NewtonianMover() {
       location = new PVector(0, 0);
     }
     void update() {
       //based on distance = initialVelocity * deltaTime + 0.5 * acceleration * time * time
       location = PVector.mult(globalAcceleration, 0.5*sq(cumulativeTime));
     }
     void display() {
       stroke(0, 255, 0);
       point(cumulativeTime, map(location.x, 0, 4000000, 0, -(height-(2*margin))));
     }
   }
   
   class EulerMover {//red line
     PVector location;
     PVector velocity;
     EulerMover() {
       location = new PVector(0, 0);
       velocity = new PVector(0, 0);
     }
     void update() {
       PVector deltaAcceleration = PVector.mult(globalAcceleration, deltaTime);
       velocity.add(deltaAcceleration);
       PVector deltaVelocity = PVector.mult(velocity, deltaTime);
       location.add(deltaVelocity);
     }
     void display() {
       stroke(255, 0, 0);
       point(cumulativeTime, map(location.x, 0, 4000000, 0, -(height-(2*margin))));
     }
   }
   
   class VerletMover {//blue line
     PVector location;
     PVector previousLocation;
     VerletMover() {
       location = new PVector(0, 0);
       previousLocation = new PVector(0, 0);//optionally comment out and replace with the alternative below
       //previousLocation = PVector.mult(globalAcceleration, -0.5 * sq(deltaTime));//this calculation is suggested in the article but seems unnecessary
     }
     void update() {
       PVector deltaLocation = PVector.sub(location, previousLocation);
       PVector deltaAcceleration = PVector.mult(globalAcceleration, sq(deltaTime));
       PVector increment = PVector.add(deltaLocation, deltaAcceleration);
       previousLocation = location.get();
       location.add(increment);
     }
     void display() {
       stroke(0, 0, 255);
       point(cumulativeTime, map(location.x, 0, 4000000, 0, -(height-(2*margin))));
     }
   }