Previously we talked about how to dampen a value to zero in a way that is independent of framerate. This time we'll look into using a spring (specifically a critically damped spring) to achieve something similar but more flexible.
Spring Overview
What is a critically damped spring? First, a spring is something that tries to maintain a certain length by exerting a force equal to the displacement from that length times a spring constant. A damped spring is one that exerts an additional force proportional to how fast its length is changing in order to counter that motion. When a spring doesn't have much (or any) damping, it oscillates. But if you get the damping just right, the spring will just decay back to its preferred length with no oscillations. This is what's called a critically damped spring.
We'll consider a theoretical spring with a point mass attached to its end at x = 0. The velocity and acceleration of the mass are x' and x'', respectively.
We start with Newton's second law of motion. The force on an object is equal to the mass of the object times the acceleration caused by the force or F = m a. As we said before, the force the spring exerts on the point mass is equal to the displacement of the spring times a constant we'll call k, in the direction opposite the displacement (so that the spring works to return the mass to x = 0).
The damping force of a spring is equal to a damping coefficient we'll call γ times the velocity of the end of the spring in the direction opposite its motion (so that the spring works to slow that motion down).
And since F = m a = m x'', then
Then we can use our knowledge about differential equations to solve for x (see the posts on linear homogeneous ODEs part 1 and part 2). Our characteristic equation is
which has roots
Now we know that if the roots are distinct and real, then our general solution will look like this
If we have one repeated root, it will look like
And if we have a complex conjugate pair of roots of the form a ± bi our solution will be
We want to avoid oscillations, so we obviously don't want the case with complex roots, since it involves sin and cos. Then, for a given spring constant k, intuitively we want to pick a damping coefficient γ as small as possible so that the spring can return to equilibrium as quickly as possible. The smallest γ can be without making the square root negative will be given by
This is the damping coefficient for a critically damped spring, and it yields the repeated root
In order to simplify much of the following, we'll create a constant
In this case of using a critically damped spring model, we don't care what the mass is relative to k, so we can just concern ourselves from here on with the constant ω. Now we have our general solution
Then if we have the initial conditions x(0) = x0 and x'(0) = x'0 we get
and
Plugging that back into our solution gives us
Figure 1 shows a few example solutions.
Figure 1: Critically damped spring movement with initial position of 1 and initial velocities of -2 (red), 0 (green) and 2 (blue)
Numerical Integration
Now we know the solution to (1), the equation that governs a spring's movement, but to use it in its current form we would need to compute c0 and c1 and then evaluate the exponential in order to update our values each frame. Depending on your situation, it might be a bit more practical to numerically integrate the position and velocity from the previous frame into an updated pair of values.
Forward Euler Integration
The simplest method is called forward Euler integration which uses the current position xn and velocity x'n to compute an acceleration and then uses that along with the time elapsed since the last update Δt to compute a new velocity x'n+1 and position xn+1. First we rewrite (1) to solve for x''
where we've substituted in the value of γ from (2). Then we use this to compute a new velocity each frame and use that in turn to compute a new position like this
This works decently well, but when you have a stiff spring (one with a large spring constant k) or if you have a frame time Δt that is very large, the change in velocity or change in position can be much too large, which can lead to your spring length oscillating or potentially blowing up and not settling at all.
Figure 2: A few iterations of a stiff spring using forward Euler integration
Figure 3: A few iterations of a non-stiff spring using forward Euler integration
Figure 2 shows an example of a few iterations using this iterator on a stiff spring with large time steps. The curves show the actual solutions from each new pair of velocity and position values so that you can clearly see how far the iterations are deviating from the intended motion. Figure 3 shows what it looks like when the spring isn't too stiff and the time steps are small enough. The numerical iteration stays fairly close to the actual solution in that case.
Backward Euler Integration
In order to deal with the potential instability we saw with forward Euler integration, let's take a look at what's known as backward Euler integration. Instead of computing acceleration from current parameters and integrating that into the velocity and position, we want our acceleration to be based off of the values at the end of the current time step instead of their current values. This means our integration looks like this
Note that this time the acceleration x''n+1 is based off of xn+1 and x'n+1 instead of xn and x'n, but these values are unknown, so we'll need to solve this as a system of equations. Substituting (5) into (6), we get
and then substituting in (7), expanding the terms and then solving for x'n+1 gives
so then one iteration of our numerical integration becomes
Figure 4: A few iterations of a stiff spring using backward Euler integration
Figure 5: A few iterations of a non-stiff spring using backward Euler integration
Figure 4 shows the first few iterations of the same stiff spring parameters from before and you can see that things are much more stable. Figure 5 shows the backward Euler integrator working on the stable spring parameters, and we can see that the backward integrator gives us less error than the forward Euler integrator.
Control
Now that we have figured out how to numerically integrate a critically damped spring such that we avoid any instability, it would be nice to have an intuitive way to set our constant ω without needing to just experiment with it until it "looks right". Let's take a hint from the post about framerate independent damping and attempt to find a similar way to specify our spring constant. In that post we specified a scalar we wanted to apply to a value every second in order to decay that value toward zero. Let's say that in this case, we want to specify the fraction of the original value remaining after a certain time period. For example, if our initial position was 5.0, we might want it to decay to 1.0 after 3 seconds, which is to say, we want to have 20 percent of the original value remaining after 3 seconds. In general we want to be able to say that we want p percent of our initial value left after s seconds.
Using the solution to the spring motion equation we got in (4) and dividing it by the initial position to get the desired remaining percent, we get
Since the remaining percent is dependent on our initial velocity, we'll simplify things by assuming a zero initial velocity. With that assumption, we get
Now we can see that solving for ω isn't easy, but we can use what's known as the Lambert W function which is the inverse of the function
So in order to use this function, we need to get our equation into a form that looks like WeW, which we can do by negating both sides and dividing by e
Since the Lambert function isn't something we can easily compute at runtime, we can set up a lookup table for the numerator of (8) and then just divide by the desired time duration s. I used Wolfram Alpha to generate my lookup table (like this).
Figure 6: Table of omega values from desired remaining percentages computed from (8) as well as empirically with code to iteratively search for the correct values
I also wrote some code to iteratively search for the numerator values based on the desired remaining fraction, and the two data sets seem to agree pretty well. Figure 6 shows the lookup table values versus the values I found empirically. The only appreciable differences are on the "stiffer" end of the spectrum. The data for figure 6 can be found here.
Conclusion
We've found that if you want to change a value smoothly toward a goal such that after s seconds the value has p fraction of the initial distance still left to travel (in other words, it has travelled 1 - p fraction of the distance toward the goal), you'll want to pre-compute ω by dividing the value you get from the above look up table by s.
Then with each update, you'll want to integrate the current velocity and current value like this
where G is the goal value (up until now we've assumed the goal was zero). It's also worth pointing out that our x values could just as easily have been vectors if you want to do more than just smoothing out scalar values.
This method of using a critically damped spring to smoothly approach a goal has been used elsewhere, but I was unable to find anyone that went into the kind of depth I wanted. I had also never seen anyone figure out a way to give the spring an intuitively controlled set of parameters, so hopefully this is useful for someone else looking for the same kind of information.
No comments:
Post a Comment