Solution to over-damped harmonic spring

Question

(A kind soul at physics.stackexchange suggested I post here as well, sorry if out of bounds.)

I'm trying to programmatically model a damped harmonic spring for use in mobile UI animations (physics mathematics isn't my background, please pardon any misconceptions).

Having derived the parameters for the general case equation, I can iteratively calculate values until I reach a suitable threshold, though because this is bound to "simple" trigonometric and $e^{x}$ functions on the CPU, the 4000-some-odd steps can cause about 0.25 seconds lag on slow devices while it calculates. I'd like to speed this up using my platform's super-optimized vector and BLAS/LAPACK variants.

The requirement for doing this is precalculating the number of steps necessary to reach my threshold value.

---

In the underdamped case, where the roots of the characteristic function of the differential equation are non-real, I can use algebraic tricks to get my values:

$$x(t) = c_{1}e^{r_{1}}\cos(i_{1}t) + c_{2}e^{r_{2}}\sin(i_{2}t)$$

(Given $r_{1}$, $i_{1}$, $r_{2}$, and $i_{2}$ are the real and irrational components of my two roots, respectively.)

Knowing that $r_{1} = r_{2}$ and $i_{1} = -i_{2}$, I can simplify to:

$$x(t) = c_{1}e^{r_{1}}\cos(i_{1}t)$$

And get my desired value of $t$ for my threshold $a$:

$$t = \arccos(a / c_{1} / e^{r_{1}}) / i_{1}$$

---

When the roots are real, the equation looks a lot simpler:

$$x(t) = c_{1}e^{r_{1}} + c_{2}e^{r_{2}}$$

However, I don't have my trig functions floating around to help me solve it (even if I did, the irrational components being 0 would cause problems, of course).

Take the concrete example on pages 3-4 of this document (my bible during this process), since they at least solve cleanly:

$$x(t) = 1.5e^{-t} - 0.5e^{-3t}$$

I know how I would solve for $t$ to get my for when $x(t) = a$ on paper, by setting $x=e^{t}$, solving, and back substituting, but I don't have that luxury here.

I can make a few assumptions: the roots and constants are all real. I'm always going to be looking for the smallest, first, positive value of $t$. Obviously, the iterative solution is the simplest for this case, but in the end that would involve more steps and therefore be slower no matter what my other optimizations would be. How, then, would I go about solving for my threshold value algorithmically in this (supposedly) simplified case?

Addendum

The underdamped solution presents an extra requirement. The motion curve will oscillate back and forth a few times across the endpoint. Therefore, "first and lowest" $t$ requirement is not necessarily true. In my current, iterative code, the threshold value is both checked against the distance from the current $x(t)$ to the endpoint, as well as to the distance from the previous $x(t)$ as well to allow for a number of oscillations. This might make a more efficient solution nearly impossible.

Answer 1

For the general equation $m \ddot{x} + b \dot{x} + k x = 0$ you can find a solution with the following steps

Updated Answer

1. From the problem coefficients calculate the characteristic values $$ \begin{aligned} \omega_n & = \sqrt{\frac{k}{m}} \\ \zeta & = \frac{b}{2 \sqrt{k m}} \end{aligned}$$
2. If $\zeta<1$ then calculate solution coefficients from the initial conditions $x(0)=x_0$ and $\dot{x}(0)=v_0$ $$ \begin{aligned} \omega & = \omega_n \sqrt{1-\zeta^2} \\ A & = x_0 \\ B & = \frac{v_0+x_0 \omega_n \zeta}{\omega} \\ r & = \omega_n \zeta \end{aligned} $$
3. The solution function is $$ x(t) = e^{-r t} \left( A \cos(\omega t) + B \sin(\omega t) \right) \tag{1}$$
4. The final time is $$ x_{tol} \approx e^{-r t_{final}} \sqrt{A^2+B^2} \\ t_{final} = \frac{ \ln{ \left( \frac{ \sqrt{A^2+B^2}}{x_{tol}} \right)} } {r} $$
5. If $\zeta>1$ then calculate the solution coefficients from the initial conditions $x(0)=x_0$ and $\dot{x}(0)=v_0$ $$ \begin{aligned} \omega & = \omega_n \sqrt{\zeta^2-1} \\ C & = \frac{1}{2} \left( x_0 + \frac{v_0 + x_0 \zeta \omega_n}{\omega} \right) \\ D & = \frac{1}{2} \left( x_0 - \frac{v_0 + x_0 \zeta \omega_n}{\omega} \right) \\ r_1 & = \left(\zeta - \sqrt{\zeta^2-1}\right) \omega_n \\ r_2 & = \left(\zeta + \sqrt{\zeta^2-1}\right) \omega_n \end{aligned} $$
6. The solution function is $$ x(t) = C e^{-r_1 t} + D e^{-r_2 t} \tag{2}$$
7. The final time is $$ \begin{aligned} x_{tol} & \approx C e^{-r_1 t_{final}} + 0 \\ t_{final} & = \frac{ \ln{\left( \frac{C}{x_{tol}} \right)}}{r_1} \end{aligned} $$ because $r_1<r_2$ and thus decaying slower.
8. For critical damping $\zeta=1$ the solution function is $$ x(t) = e^{-\omega_n t} \left(x_0 + (v_0+x_0 \omega_n) t \right) \tag{3}$$
9. The final time is $$ \begin{aligned} \ln{\left(x_{tol}\right)} & \approx {\rm Taylor}\left( \ln{\left( (v_0 + x_0 \omega_n) t + x_0 \right)} -\omega_n t, t=\frac{2\pi}{\omega_n},{\rm order}:=1 \right) \\ X & = \frac{2\pi v_0}{\omega_n} + x_0 (2\pi+1) \\ t_{final} &= \frac{ 2 \pi \left( X \ln \left( \frac{X}{x_{tol}} \right) -X+x_0 \right)}{\omega_n \left( X ( 2\pi -1)+x_0 \right)} \end{aligned}$$

Answer 2

As your document says, when you have two real roots and the system is overdamped, the one with the longer time constant (inverse of the root) controls how close you are to zero. So if you are trying to see how long it takes to get close to zero (if $a$ is much smaller than the original displacement) you can ignore the one that falls fastest. In your example, that is the $e^{-3t}$ term. You can then solve $a=1.5e^{-t}$ to get $t=-\log \frac a{1.5}$. As an example, if $a=0.1$, we get $t \approx 2.7$ The other term then contributes $-0.00015$, quite negligible.

If they are very close, you can use Newton's method or fixed point iteration. Say you were solving $0.1=1.5e^{-t}-0.5e^{-1.01t}$ You would be looking for a root of $f(x)=-0.1+1.5e^{-t}-0.5e^{-1.01t}$, which gives $f'(x)=-1.5e^{-t}+0.505e^{-1.01t}$ For a starting value, ignore the difference of the exponentials and take $x_0=-\log \frac {0.1}{1.5-0.5}\approx2.3$ and $x_{i+1}=x_i-\frac {f(x_i)}{f'(x_i)}$ which will converge quickly.