Game Math: Precise Control over Numeric Springing

This post is part of my Game Math Series.

Source files are on GitHub

Check out this post if you want to see more visual examples of numeric springing.

Numeric springing is a very powerful tool for procedural animation. You specify the initial value, initial velocity, target value, and some spring-related parameters; the result is a smooth springing effect. You can apply this technique to all sorts of numeric properties, some common ones being position, rotation, and scale of an object.

spring

The Common But Designer-Unfriendly Way

Here is a common, yet designer-unfriendly, implementation of numeric springing. Let x_i and v_i denote the numeric value and velocity, respectively, at frame i; and let k, d, h, x_t, denote the spring stiffness, damping factor, time step between frames, and target value, respectively.

     \begin{flalign*}   v_{i+1} &= d v_i + h k (x_t - x_i) \\   x_{i+1} &= x_i + h v_{i+1} \\ \end{flalign*}

Let’s dissect this implementation.

The difference between the target value and the current value (x_t - x_i) contributes to the velocity change after multiplied by the spring stiffness k. The larger the difference and spring stiffness, the larger the change in velocity.

The new velocity v_{i+1} inherits the value from the old velocity v_i scaled by the damping factor d, which should be a value between zero and one. If the damping factor is equal to one, then no scaling happens, and the new velocity completely inherits the value from the old velocity. If the damping factor is a fraction of one, then the magnitude of oscillation will decrease gradually.

That seems all reasonable and good, but just what exactly are the values for k and d do we need to achieve a specific springing effect, say an oscillation at 5Hz frequency and the oscillation magnitude decreases 90% every one second? Many people just decide to bite the bullet and spend a large amount of time hand-tweaking the values for k and d, until the result is acceptably close to the desired effect.

We Can Be Exact

Now let’s take a look at the differential equation for a damped spring system centered around x = x_t.

    \[ \frac{\mathrm{d}^2 x}{\mathrm{d} t^2}} + 2 \zeta \omega \frac{\mathrm{d} x}{\mathrm{d} t}} + \omega^2 (x - x_t) = 0 \]

I won’t go into details as to how this equation is derived. All you really need to know is that \omega (omega) and \zeta (zeta) govern how the system behaves and what they mean.

\omega is the angular frequency of the oscillation. An angular frequency of 2\pi (radians per second) means the oscillation completes one full period over one second, i.e. 1Hz.

\zeta is the damping ratio. A damping ratio of zero means there is no damping at all, and the oscillation just continues indefinitely. A damping ratio between zero and one means the spring system is underdamped; oscillation happens, and the magnitude of oscillation decreases exponentially over time. A damping ratio of 1 signifies a critically damped system, where this is the point the system stops showing oscillation, but only converging to the target value exponentially. Any damping ratio above 1 means the system is overdamped, and the effect of springing becomes more draggy as the damping ratio increases.

Here’s a figure I borrowed form Erin Catto‘s presentation on soft constraints, to show you the comparison of undamped, underdamped, critically damped, and overdamped systems.

damping ratios

Below are the equations for simulating a damped spring system using the implicit Euler method.

     \begin{flalign*}   x_{i+1} &= x_i + h v_{i+1} \\   v_{i+1} &= v_i - h (2 \zeta \omega v_{i+1}) - h \omega^2 (x_{i + 1} - x_t) \\ \end{flalign*}

Solving x_{i+1} and v_{i+1} in terms of x_i and v_i using Cramer’s rule, we get:

     \begin{flalign*}   x_{i+1} &= \frac{\Delta_x}{\Delta} \\   v_{i+1} &= \frac{\Delta_v}{\Delta},  \\ \end{flalign*}

where:

     \begin{flalign*}   \Delta &= (1 + 2 h \zeta \omega) + h^2 \omega^2 \\   \Delta_x &= (1 + 2 h \zeta \omega) x_i + h v_i + h^2 \omega^2 x_t \\   \Delta_v &= v_i + h \omega^2 (x_t - x_i) \\ \end{flalign*}

Below is a sample implementation in C++. The variables x and v are initialized once and then passed into the function by reference every frame, where the function keeps updating their values every time it’s called.

/*
  x     - value             (input/output)
  v     - velocity          (input/output)
  xt    - target value      (input)
  zeta  - damping ratio     (input)
  omega - angular frequency (input)
  h     - time step         (input)
*/
void Spring
(
  float &x, float &v, float xt, 
  float zeta, float omega, float h
)
{
  const float f = 1.0f + 2.0f * h * zeta * omega;
  const float oo = omega * omega;
  const float hoo = h * oo;
  const float hhoo = h * hoo;
  const float detInv = 1.0f / (f + hhoo);
  const float detX = f * x + h * v + hhoo * xt;
  const float detV = v + hoo * (xt - x);
  x = detX * detInv;
  v = detV * detInv;
}

Designer-Friendly Parameters

Now we have our formula and implementation all worked out, what parameters should we expose to the designers from our spring system so that they can easily tweak the system?

Inspired by the aforementioned presentation by Erin Catto, I propose exposing the oscillation frequency f in Hz, and the fraction of oscillation magnitude reduced p_d over a specific duration t_d due to damping.

Mapping f to \omega is easy, you just multiply it by 2 \pi:

    \[ \omega = 2 \pi f \]

Mapping p_d and t_d to \zeta is a little bit more tricky.

First you need to understand that the oscillation magnitude decreases exponentially with this curve:

    \[ y(t) = e^{-\zeta \omega t} \]

After plugging in p_d and t_d, we get:

    \[ y(t_d) = e^{-\zeta \omega t_d} = p_d \]

So we can solve for \zeta:

    \[ \zeta = \frac{\ln (p_d)}{-\omega t_d} \]

For example, if we want the oscillation to reduce 90% (p_d = 0.1) every 0.5 second (t_d = 0.5) with oscillation frequency 2Hz (\omega = 4 \pi), then we get:

    \[ \zeta = \frac{\ln (0.1)}{-4 \pi \cdot 0.5} = 0.367 \]

Conclusion

Numeric springing is a very powerful tool for procedural animation, and now you know how to precisely control it. Given the desired oscillation frequency and percentage of oscillation magnitude reduction over a specific duration, you can now compute the exact angular frequency \omega and damping ratio \zeta necessary to create the desired springing effect.

By the way, let’s look at the animated example shown above again:

spring

This animation was created with \omega = 8 \pi (4 full oscillation periods per second) and \zeta = 0.23 (oscillation magnitude reduces 99.7% every second).

About Allen Chou

Physics / Graphics / Procedural Animation / Visuals
This entry was posted in C/C++, Gamedev, Math. Bookmark the permalink.

9 Responses to Game Math: Precise Control over Numeric Springing

  1. Chris Wade says:

    Hey Allen, thanks for writing this article. It’s leads helping the tweens in my game. One suggestion: the code snippet having abbreviated argument names like x, v, xt, etc made the article and especially the spring function much harder to parse. I think it would help new readers understand quickly if those function had full names.

    • Allen Chou says:

      Thanks for the feedback. I’m glad you find it helpful.
      And yes, I did consider using variable names that are more descriptive; however, I decided that it’s more important to have the code nicely formatted without line-wrapping in the blog post, so I kept the variable names short. I think the comments above the function explaining the variables should be enough.

  2. Pingback: Precise Control over Numeric Springing |

  3. davidtheory says:

    Hi Allen, thanks for the post! I just wanted to add one cent to this. If the oscillator has a maximum amplitude of M at start t=0 (in the post M=1), then we just need to multiply Pd by M.

    • Allen Chou says:

      Can you explain? I thought Pd already specifies the percentage of reduction (regardless of M), so M is reduced to Pd * M after duration td. Or am I missing something?

      • davidtheory says:

        What I meant is that in the post |x(0)| = 1, but what if I want |x(0)| = M, then |x(td)| = Pd*M. However, now that I think more about it, it is better to leave x(t) to be in [-1,1] as it is in your post, and use its value to whatever amplitude we need outside (just like we do with easing functions). Sorry for the confusion :/

        • Allen Chou says:

          Actually, it’s only the figure that says |x(0)| = 1. The equations I showed don’t require |x(0)| = 1. Perhaps I should have clarified that.

          • davidtheory says:

            wait, but we have explicitly x(t) = exp(-zwt), and from here we have x(t=0) = exp(0) = 1 , or am I missing something?

          • Allen Chou says:

            Yeah. That’s a typo. I meant to write y(t) = exp(-zwt), representing the oscillation reduction percentage (fixed now).

Leave a Reply