6
\$\begingroup\$

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.
enter image description here
(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

\$\endgroup\$
  • \$\begingroup\$ What level maths do you know? \$\endgroup\$ – MickLH Aug 16 '15 at 15:30
  • \$\begingroup\$ @MickLH, i'm a University student at Software Engineering major, i'm not that good so if its not something super advance i should be able to understand it. \$\endgroup\$ – ali ahmadi Aug 16 '15 at 15:43
6
\$\begingroup\$

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:

enter image description here

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:

enter image description here

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:

enter image description here

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

enter image description here

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:

enter image description here

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:

enter image description here

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, [ ..... ]]);
\$\endgroup\$
  • \$\begingroup\$ i will try to implant this by tomorrow , hope it works, for now it looks like something that can work :) so Thank You. \$\endgroup\$ – ali ahmadi Aug 16 '15 at 18:37
  • \$\begingroup\$ i have a few question 1- where should i use friction [n+1] ? 2- in friction equation , you have used friction ! when u say pow(e, (-decayTime * fiction) / 2) 3- for me to get the current position , i should get y(t) right ? i mean i give t , the passedTime , right ? \$\endgroup\$ – ali ahmadi Aug 17 '15 at 13:05
  • \$\begingroup\$ if possible, can you copy ur code here ? i cant get the friction right, if possible i'd like to analyise your code to see how its implanted, thx \$\endgroup\$ – ali ahmadi Aug 17 '15 at 13:53
  • \$\begingroup\$ please check my "Edit" in the main post, i still have some problem \$\endgroup\$ – ali ahmadi Aug 18 '15 at 14:16
  • \$\begingroup\$ is it possible that because you are using maxima , somehow it fixes the divide on zeros and ignores them ? and now that i'm implanting them on c++ its causing problem ?, i know you have spend quite some time on your answer and that is exactly why i m asking you to help me one more time so i can get some results, please check my main post, i just cant get the friction right, Thank you for your help. \$\endgroup\$ – ali ahmadi Aug 20 '15 at 9:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.