physics - How can I fix my velocity damping to work with any delta/frame time?

I am decreasing my velocity by 50% every second using a guide I found online. I am using the code here and it gives the right result but only for very small values for dt:

float dt, velocity, position;

dt = 0.002f;

velocity = 1.0f;
position = 0.0f;

for (float t = 0.0f; t <= 1.0f; t += dt) {
    position += velocity * dt;
    velocity = lerp(velocity, 0.0f, 1.0f - std::pow(0.5f, dt));
}

printf("dt = %f, p = %f, v = %f\n", dt, position, velocity);


dt = 0.09f;
velocity = 1.0f;
position = 0.0f;

for (float t = 0.0f; t <= 1.0f; t += dt) {
    position += velocity * dt;
    velocity = lerp(velocity, 0.0f, 1.0f - std::pow(0.5f, dt));
}

printf("dt = %f, p = %f, v = %f\n", dt, position, velocity);


return 0;

But I am getting roughly the right velocity BUT two different final positions back when I use two different delta times (not due to floating point issues):

dt = 0.002000, p = 0.722848, v = 0.499308
dt = 0.090000, p = 0.784219, v = 0.473029

How can I fix this difference inside of my update loop?

How do I account for the extra acceleration between timesteps?

Answer

Although you're correcting your velocity, your position update:

position += velocity * dt;

...is still treating velocity as though it's a constant over the frame duration. So, the smaller your timestep, the more accurate this simplification is, but with longer timesteps, it gives too much weight to the high velocity early in the frame, without accounting for how it diminishes later in the frame.

To get this right, we'll need to integrate velocity over the frame duration with some calculus.

If your velocity is equal to...

$$\begin{align}\vec v_{\Delta T} &= \vec v_0 \cdot r^{\Delta T} \\&=\vec v_0 \cdot e^{\Delta T \ln r}\end{align}$$

Then your position is the integral of this value over your timestep:

$$\begin{align} \vec p_{\Delta T} &= \vec p_0 + \int_0^{\Delta T} {\vec v_0 \cdot e^{t\ln r}\ dt} \\&= \vec p_0 + \vec v_0 \cdot \left(\frac {e^{t \ln r}} {\ln r} \right) \Bigg|_0^{\Delta T} \\&= \vec p_0 + \vec v_0 \cdot \left(\frac {r^{\Delta T} - 1} {\ln r} \right) \end{align}$$

In code:

position += velocity * (std::pow(r, dt) - 1.0f)/std::log(r);
velocity *= std::pow(r, dt);