How To Program a fully controllable spring (Damped) motion

Question

i'm working on a project that needs realistic 2D Animation, none linear motion you might say, i have tried many things such as Hermite Curve and a few other motions and in the end Spring Motion.
because my spring is based on a simple distance, velocity, friction base i can not control the time that spring finishes and the number of times it pases through it's target, after some study i found out that i need to program damping motion and i found an article
http://www.ncsu.edu/per/Articles/MarchewkaAbbott&Beichner.pdf
which uses Total energy to solve the problem, but even after reading it many times i couldn't program it !
which is the reason i'm here , please i really need help
my goals : programming a fully controllable damping motion that is
1- has a progress based on an interval between 0 and 1 (meaning whatever the motion is, you can play it like a movie, you give the percent and it returns the position)
2- you can determine the number of times that object passes through the target before it settles down here is the final function that i have in mind :

 /*Returns the position object must be at the specified percent */  
 springAtPercent(startPosition, targetPosition, cycleNum /*the Number Of Times That Object Passes Through The Target*/, percent /*between 0.0 and 1.0*/);

i think the function makes it prettly clear what i need, i have searched a lot and found lots of result ! but none of them is good enough or has enough maneuver capability to achieve complicated animations that is why i need something like this.
in the end here is an image i made for a question in StackExchange physics section which resulted in the article i mentioned, hope it helps to further clarify my question.

(i'm programming in pure C++, i also use ogre, in the end programming language doesn't matter that much, so if you don't know c++ feel free to use any language you like)
Thanks.

Edit 1:
i have tried to implant MickLH answer
but i have a few problem!
i can not get the friction right ! i know what Newton's Method is and how it works but how many times do i need to calculate the new friction ?
and when you say friction[0] = 0 do you mean the x(n+1) = x + f(x) / f'(x) with n = 0 , should have an x = 0 ? because i tried that and the result of friction is very small! (as you can see the code)
i think the stiffness formula is right ! because i tried with test variable(which is 7.164850514373732 from your code) and with that constant i can get a right stiffness, which means the only problem i have is with the friction ! can you explain what should happen so that i get a right friction ? (note i can not use any math library because i'm using an company engine), thx

the lines i have written for friction and stiffness :

// as you can see the very first friction is calculated using friction = 0 which means e ^ 0 = 1 
ATVector2 friction = csc(Ogre::Math::PI * cycleNum) * (((2 * Ogre::Math::PI * cycleNum * Ogre::Math::Cos(Ogre::Math::PI * cycleNum)) / decayTime) - ((2 * Ogre::Math::PI * cycleNum * 1 * pull) / (decayTime * stopped)));
    ATVector2 aftePow =  ATVector2(pow(e, -(decayTime * friction.x) / 2), pow(e, -(decayTime * friction.y) / 2));
    ATVector2 frictionplus1Part1 = Ogre::Math::PI * cycleNum * (2 * Ogre::Math::Cos(Ogre::Math::PI *  cycleNum) * aftePow * stopped - decayTime * friction * pull - 2 * pull);
    ATVector2 frictionplus1Part2 = decayTime * (Ogre::Math::Sin(Ogre::Math::PI * cycleNum) * aftePow * stopped - (Ogre::Math::PI * cycleNum * pull));
    ATVector2 frictionplus1 = frictionplus1Part1 / frictionplus1Part2;
    ATVector2 test = ATVector2(1.842068074395237, 1.842068074395237);
    ATVector2 stiffness = ((decayTime * decayTime) * (test * test) + 4 * (Ogre::Math::PI * Ogre::Math::PI) * (cycleNum * cycleNum)) 
    / (4 * decayTime * decayTime);

there is also csc which is 1 / sin(pi *4) which is divide on 0 , which if we fix, can be a large number about 672, is this okay as well ?
here is another test which even after 100 time calculating friction i can only get 1.#INF0000 which is invalid!
here is the test :

//here, the initial value for friction is set to zero
mFriction = ATVector2::ZERO;
ATVector2 aftePow = ATVector2::ZERO;
/*
    decayTime = 5;
    cycleNum = 4;
    stopped = 0.01;
    pull = ATVector(1, 1);
*/

for (int i = 0; i < 100; i++)
{

    // here i calculate e^ ... to make things simpler
    aftePow =  ATVector2(pow(e, -(decayTime * mFriction.x) / 2), pow(e, -(decayTime * mFriction.y) / 2));
    mFriction = csc(Ogre::Math::PI * cycleNum) * (((2 * Ogre::Math::PI * cycleNum * Ogre::Math::Cos(Ogre::Math::PI * cycleNum)) / decayTime) - ((2 * Ogre::Math::PI * cycleNum * aftePow * pull) / (decayTime * stopped)));
}

//after the loop mFriction = ATVector(1.#INF0000, 1.#INF0000)

Edit 2 :
for now the friction is not going to be an issue! because i found an even bigger problem!
in MickLH answer, the friction is calculated on the onUpdate() and it uses the previous friction and that makes it useless inside a game with delta time on each update, why?
for example when decayTime is 5 it means if the animation is played for example 10 time in every second, it will play a smooth and good spring! but if some kind of lag occurs that means the friction was not updated correctly! and the animation will be played on decayTime + lagTime !
it also take away control of the element for example :
you can not play an specific part of the animation, for example in a 5 sec animation you only need the last 2 second ! using MickLH answer you can only do that if you write some kind of loop and do the calculation 10 * 3 = 30 time ! so the friction is on the right amount !
in the end The Answer does not meet my need for getting the position on a specific percent ! , which is also something i mentioned at the beginning of my question
so once again if possible, anyone, please i need a function that can return a position based on Percent, based on DeltaTime, based on the time that passes, something accurate , i say this once again, it must become something like
this :

springAtPercent(startPosition, targetPosition, cycleNum /*the Number Of Times That Object Passes Through The Target*/, percent /*between 0.0 and 1.0*/);  

MickLH is it possible to write your answer based on DeltaTime, percent ?
Thx

Answer 1

I'm using Hooke's Law here as the definition of a spring. (\$F = -k * X\$)

Given the derivatives of position and velocity, are velocity and force respectively, we can construct a differential equation for the stretching of the spring.

$$ \frac{d^2}{dt^2}y(t)=-\textit{Stiffness }y(t) - \mathit{Friction} \bigg{(}\frac{d}{dt}y(t)\bigg{)}$$

Which is just a damped harmonic oscillator, and since we already know that only the under-damped case need analysis, we can obtain a nice solution:

$$ \small{ y(t) = e^{-\frac{\textit{Friction }t}{2}}\bigg{(}\%k_1 \text{ sin}\big(\frac{t\sqrt{4\textit{ Stiffness} - Friction^2}}{2}\big) + \%k_2 \text{ cos} \big( \frac{t\sqrt{4 \textit{ Stiffness} - \textit{Friction}^2}}{2} \big) \bigg)} $$

We can eliminate the unknowns by knowing initial position \$ y(0) = pull \$ and velocity \$ \frac{d}{dt}y(0)=0 \$

$$ \small{ y(t) = \frac { pull \sqrt{4\:\textit{Stiffness} - \textit{Friction}^2}} { \scriptsize{ \sqrt{4\:\text{Stiffness} - \textit{Friction}^2}\:e^\frac{Friction\:t}{2} \text{cos}\big(\frac{\sqrt{4\:\text{Stiffness} - \text{Friction}^2} t}{2}\big) - Friction\;e^\frac{Friction\;t}{2} \text{sin} \big( \frac{\sqrt{4\:\text{Stiffness} - \text{Friction}^2} t}{2} \big) }}}$$

Ugly. But we have a closed form solution in time :)

Before we can solve for the parameters you're interested in, we have to address one small ambiguity: The damped harmonic oscillator never stops, only decays. I will use a threshold where we consider motion "stopped", and solve for the peak which attains this amplitude.

Now, from the solution above, I have obtained that the set of peaks are generated by:

$$t=\frac{2\pi\mathbb{Z}}{\sqrt{4\text{ Stiffness} - \text{Friction}^2}}$$

Which conveniently brings the n-th peak into the relationship.

Since you've specified both the decay time, and the number of crossings, we can derive a relationship between the two constants we are trying to compute:

$$\text{Stiffness}=\frac{\text{DecayTime}^2\text{Friction}^2+4\pi^2\text{CycleNum}^2}{4\text{ DecayTime}^2}$$

And now by substituting the time of the peak into the amplitude function, we can express the magnitude of the peak after a target number of equilibrium crossings:

Setting this equal to our "Stopped" threshold, we can then solve for the desired Friction coefficient:

I have obtained an implicit solution, and although it looks like it could maybe be solved by the W function in closed form, we would then need to implement the W function in C++, so I have not investigated that. Instead, I have taken the approach of reformulating the problem as root-finding, and numerically approximating the soluiton using Newton's Method. This strategy has yielded the following fixed point iteration:

I have not analysed the convergence properties of this formula with respect to the initial guess. With that said, seeding with Friction[0] = 0 appears to place the attractor into the correct basin.

As a quick test, I've plugged in these values:

CycleNum = 4
DecayTime = 5
Stopped = 0.01
Pull = 1

And received the results:

Stiffness = 7.164850514373732
Friction = 1.842068074395237

Which I have then plugged into the most basic Euler integrator, with a step size of 1/10, and it produced this graph:

Which appears to be right on target, crossing 4 times to achieve a peak value of 0.01, after 5 seconds of simulated time.

---

EDIT

Here's my prototype code, I let maxima do newton's method:

params:[DecayTime = 5.0, Stopped = 0.01, Pull = 1.0, CycleNum = 4.0]$

load("mnewton")$
Friction = csc(%pi*CycleNum)*(2*%pi*CycleNum*cos(%pi*CycleNum)/DecayTime
                              -2*%pi*CycleNum*%e^-(DecayTime*Friction/2)*Pull
                               /(DecayTime*Stopped))$
subst(params, %)$
mnewton(%, Friction, 0);

[[Friction=1.842068074395236]]

I also continued with maxima while completing stiffness:

subst(flatten([%,params]), Stiffness = (DecayTime^2*Friction^2+4*%pi^2*CycleNum^2)/(4*DecayTime^2)), numer;

Stiffness=7.164850514373732

So I plugged these constants into a quick little C++ program:

#include <stdio.h>

class SpringAndMassEuler {
public:
    float Stiffness, Friction;
    float MassOffset, MassVelocity;

    SpringAndMassEuler() {
        Stiffness = 7.164850514373732;
        Friction = 1.842068074395237;
        MassOffset = 1;
        MassVelocity = 0;
    }
    
    float Update(float dt) {
        MassVelocity -= dt * (Stiffness * MassOffset + Friction * MassVelocity);
        MassOffset += dt * MassVelocity;
        return MassOffset;
    }
};

int main() {
    float t = 0;
    
    SpringAndMassEuler sim;
    printf("[%f, %f],\n", t, sim.MassOffset);
    
    // Choosing this timestep to play it safe with lame Euler
    float dt = 0.1;
    
    int c = 50;
    while (c--) {
        sim.Update(dt);
        t+=dt;
        printf("[%f, %f],\n", t, sim.MassOffset);
    }
}

Which produced the data pairs that I copy pasted back into maxima:

plot2d([discrete, [ ..... ]]);