1
\$\begingroup\$

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.

\$\endgroup\$

2 Answers 2

1
\$\begingroup\$

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.

\$\endgroup\$
2
  • \$\begingroup\$ Using the code above to compute s at t=1 where acceleration is 1 and timestep is 0.1 I get s = 0.55. If I lower the timestep to 0.01, then the result gets worse: s = 0.50 when the correct value, according to the displacement formula is 1. I am lowering the timestep and also geting worse results? That goes against what you said. Why is that? \$\endgroup\$ Jul 8, 2019 at 23:29
  • \$\begingroup\$ oh I'm tripping, with timestamp=0.01, at t=1 displacement is only 0.5 and t=2 is 2.01, so you're right, the smaller the timestamp the more accurate it is. Awesome, got it. Thanks!! \$\endgroup\$ Jul 9, 2019 at 0:36
1
\$\begingroup\$

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
\$\endgroup\$
7
  • \$\begingroup\$ I think your example is not correct. s(20 steps): 2.100000 should have been 2.0 according to the displacement formula: s = u*t+0.5f*a*pow(t, 2.0f); where t=2 = s = 0*2 + 0.5f * 1.0 * pow(2, 2.0f) = 2.0; My question is exactly about that 0.1 at t = 2. Why is it 2.1 and not 2.0 as the formula correctly outputs? Notice there are no floating point errors at t=2; \$\endgroup\$ Jul 8, 2019 at 23:09
  • \$\begingroup\$ @DaniloSouzaMorães look closely at the code, first we calculate exact solution and output it first, then we increment steps three times which increases accuracy of the integration, and finally I show that in case of single precision floating point type, at extreme amount of steps happens huge error accumulation, which is not the case for double precision. \$\endgroup\$
    – Ocelot
    Jul 9, 2019 at 17:03
  • \$\begingroup\$ My example clearly stated timestep = 0.1 as being one of the values I couldn't explain. You answer doesn't address that one. Single precision isn't the culprit for those values. It is a valuable answer nonetheless, hence the upvote. \$\endgroup\$ Jul 10, 2019 at 21:52
  • \$\begingroup\$ @DaniloSouzaMorães my answer does address that one, it's when we have 20 steps. And the explanation for that is literally the whole code, to show you that the error becomes smaller with the increasing amout of steps. But I guess I should have just linked a few wikipedia articles instead and called it a day. \$\endgroup\$
    – Ocelot
    Jul 11, 2019 at 9:01
  • \$\begingroup\$ You dont get it. At 20 steps, both double and float precision result in 2.1 when the correct answer according to the displacement formula is 2. That 0.1 in error is due to the Euler method as Alex pointed out, not due to floating point precision. \$\endgroup\$ Jul 11, 2019 at 18:00

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .