Calculate position given velocity and acceleration at a given timestep

Question

I am following a Pixar course at Khan Academy and I came across a simulation of a double spring with a given timestep. I understand how it works and why it behaves the way it does, but I want to know if the math is actually accurate. What I'm interested in is the following lines:

 // Mass 1 velocity
 mass1VelocityY = mass1VelocityY + mass1AccelerationY * timeStep;
 mass1VelocityX = mass1VelocityX + mass1AccelerationX * timeStep;

 // Mass 1 position
 mass1PositionY = mass1PositionY + mass1VelocityY * timeStep;
 mass1PositionX = mass1PositionX + mass1VelocityX * timeStep;

So we first calculate the velocity, adding a fraction of the calculated acceleration to it, and then we calculate the position, adding a fraction of the calculated velocity to it. So we are building the final position at time T by summing T/timestep intermediate positions.

Given the formula for displacement s=ut+0.5at^2 where:

s = displacement u = initial velocity a = acceleration t = time

If I plug in the values: s = 0, u = 0, a = 1, t = 2, I will get s = 2. That means that if I start at position zero and I build up the velocity with an acceleration of 1m/s/s I will end up at position 2, with a velocity of 2m/s.

Now if I try to follow the same logic but break that result into 20 steps(2/0.1 - t = 2 timeStep = 0.1), and sum the results of all these intermediate steps as is being done in that code, I will get a different result: s=0.55 at t=1 and s=2.1 at t=2.

My initial intuition is that because in the code we end up multiplying timeStep by acceleration twice, it becomes exponentiation and so the progession isn't linear anymore and so smaller steps will get smaller values in the beginning. So I have 3 questions:

1) Did I understand what the code is doing correctly?

2) Is what the code is doing the most correct/accurate way to calculate displacement?

3) I really want to grasp these concepts so if you have any other advice or know of something that would be helpful for me to learn, please point me to it.

Answer 1

This is the “explicit” or “forward” Euler method for numerically integrating your differential equation that models a spring.

https://en.wikipedia.org/wiki/Euler_method

It is unconditionally unstable, meaning that no matter what time step you use, the error with respect to the true solution to the equations can compound without bound as time goes on.

The smaller the time step, however, the slower the error accumulates. That’s what you’re seeing here.

Another approach is the “implicit” or “backward” Euler method, which is unconditionally stable. This computes a solution whose error does not compound with abandon like that of the forward approach, even if the time steps are large.

https://en.m.wikipedia.org/wiki/Backward_Euler_method

This is a deep topic of interest to all sorts of people, especially those with practical interest in fast and correct solutions to various differential equations. Those people have built up a lot of theory about it over the years.

Answer 2

No, that's because of the nature of mathematical integration and floating point precision. If you break it into more steps you will get more accurate results because of the better integration, but if you make too much steps, you will actually start to get more errors because of the insufficient floating point precision.

#include <cmath>
#include <cstdio>

//uncomment to test double precision
//#define DOUBLE_PRECISION

#ifdef DOUBLE_PRECISION
typedef double fp_type;
#else
typedef float fp_type;
#endif

void run_test(int steps, fp_type u, fp_type a, fp_type t)
{
    fp_type timestep = t / steps;
    fp_type mass1Velocity = 0.0f;
    fp_type s = 0.0f;
    for(int i = 0; i < steps; i++)
    {
        mass1Velocity += a * timestep;
        s += mass1Velocity * timestep;
    }
    printf("s(%i steps): %f\n", steps, s);
}

int main()
{
    //initial velocity
    fp_type u = 0.0f;
    //acceleration
    fp_type a = 1.0f;
    //time
    fp_type t = 2.0f;
    //displacement
    fp_type s = u*t+0.5f*a*pow(t, 2.0f);
    printf("s(exact solution): %f\n", s);

    //test 2e(i) steps
    for(int i = 0; i < 7; i++)
        run_test(2*pow(10, i), u, a, t);

    return 0;
}

Output (single precision):

s(exact solution): 2.000000
s(2 steps): 3.000000
s(20 steps): 2.100000
s(200 steps): 2.009999
s(2000 steps): 2.001011
s(20000 steps): 2.000240
s(200000 steps): 2.002040
s(2000000 steps): 1.987711

Output (double precision):

s(exact solution): 2.000000
s(2 steps): 3.000000
s(20 steps): 2.100000
s(200 steps): 2.010000
s(2000 steps): 2.001000
s(20000 steps): 2.000100
s(200000 steps): 2.000010
s(2000000 steps): 2.000001