Pendulum movement is half done

Question

I'm studying to implement a pendulum motion simulation.

void SinglePendulum::setPendulum(double dt) {
    // Apply gravity
    auto acc = acceleration;
    Vec3<double> gravity(0.0, -9.8, 0.0);   // x = 0.0, y = -9.8, z = 0.0
    acc += gravity * inverseMass;           // get accelation
    velocity += acc * dt;                   // get velocity
    velocity *= pow(dampingCoeff, dt);      // apply damping
    variablePosition += velocity * dt;      // gravity-affected weight

    // Vector from reference point to end point after being affected by gravity
    variableVector = variablePosition - pivotPosition;        // a Position - b Position : b->a vector

    // The angle between the reference vector and the current vector
    double angle;
    angle = standardVector.Dot(variableVector);                             // dot product the reference vector and the variable vector
    double standVecLength = standardVector.Length();                        // magnitude of reference vector
    double variaVecLength = variableVector.Length();                        // magnitude of variable vector
    double finalAngle = acos(angle / (standVecLength * variaVecLength));    // get angle

    // Apply position with pendulum motion
    changePos.Set(sin(finalAngle) * length, cos(finalAngle) * length, 0.0); // apply x, y, z value
}

As far as I know, It is understood that after gravity, the pendulum motion can be implemented using the angle between the reference straight vector and the changed vector.

So I apply gravity first, I then calculated the vector from the reference point to the point affected by gravity.

And I did dot product with the reference vector and the vector of the changed position.

And I calculated the magnitude of each vector, and I calculated the angle using magnitude value and theta value.

And the angle is sin on the x-axis and cos on the y-axis, which is the pendulum motion.

However, the following results were found.

I have no idea which part I missed. Is there anything else I should do?

I looked at this document and tried to modify the code, but it didn't work properly either.

Edited code

void SinglePendulum::setPendulum(double dt) {
    // The angle between the reference vector and the current vector
    double theta;
    theta = standardVector.Dot(variableVector);                             // dot product the reference vector and the variable vector
    double standVecLength = standardVector.Length();                        // magnitude of reference vector
    double variaVecLength = variableVector.Length();                        // magnitude of variable vector
    double angle = acos(theta / (standVecLength * variaVecLength));         // get angle

    Vec3<double> gravity(0.0, -9.8, 0.0);   // x = 0.0, y = -9.8, z = 0.0

    Vec3<double> forceX;
    Vec3<double> forceY;
    forceX.SetX(-9.8 * sin(angle));     // x = -9.8 * sin(angle), y = 0.0, z = 0.0
    forceY.SetY(-9.8 * cos(angle));     // x = 0.0, y = -9.8 * cos(angle), z = 0.0

    // Apply gravity
    auto acc = acceleration;
    acc += gravity * inverseMass;           // get accelation
    acc += forceX;
    acc += forceY;
    velocity += acc * dt;                   // get velocity
    velocity *= pow(dampingCoeff, dt);      // apply damping
    variablePosition += velocity * dt;      // gravity-affected weight

    // Vector from reference point to end point after being affected by gravity
    variableVector = variablePosition - pivotPosition;      // a Position - b Position : b->a vector
}

Answer 1

From mine point of view you are doing the sim wrongly ...

• acc,vel should be also vec3 (not seen in your code so not sure if it OK or not)
• no need for goniometrics vector math is enough (acos is not sign safe)
• your dampening is not how friction works its magnitue depend on actual velocity (this is the major reason of the problem you have)

This is how I see it should be (C++, GLSL like vector math):

//---------------------------------------------------------------------------
vec3 p0;            // fixed point of pendulum
float l;            // pendulum length (from p0 to pos)
vec3 pos,vel;       // ball state
float r;            // ball radius
//---------------------------------------------------------------------------
void pendulum_update(double dt)
    {
    vec3 acc,dir;
    float v;
    const float k=0.0001,g=9.81;
    dt*=10.0;                                   // time multiplier for simulation speed ...
    // compute driving force/acceleration
    acc=vec3(0.0,+g,0.0)        // gravity
       -vel*(k*length(vel));    // air friction
    // remove radial part of acc
    dir=pos-p0; dir/=l;         // unit direction vector of pendulum
    acc-=dir*dot(dir,acc);
    // Newton/d'Lambert simulation
    vel+=acc*dt;
    pos+=vel*dt;
    // repair pos so its fixed distance to p0
    dir=pos-p0; dir/=length(dir);               // unit direction vector of pendulum
    pos=p0+l*dir;
    // repair vel so its perpendicular to pendulum
    v=length(vel);
    if (v>1e-6)
        {
        vel=cross(vel,dir);
        vel=cross(dir,vel);
        vel*=v/length(vel);
        }
    }
//---------------------------------------------------------------------------
void pendulum_init(float xs,float ys)               // init pendulm to screen with xs,ys dimensions
    {
    Randomize();
    r=0.025*xs;
    l=0.25*xs;
    p0 =vec3(0.5*xs,0.5*ys,0.0);
    pos=p0+vec3(+l,0.0,0.0);
    vel=vec3(0.0,0.0,0.0);
    }
//---------------------------------------------------------------------------

I just slightly modified my Can't flip direction of ball without messing up gravity example to match yours.

Here preview:

[Edit1] double pendulum

Not sure if physically correct but I see it like this:

//---------------------------------------------------------------------------
struct _ball            // mass point
    {
    int ix;             // (-2): fixed, (-1): free, (0,1,2...): joined to ball[ix] using l sized solid bond, no loops allowed!!!
    float r,m,l;        // ball radius [m], ball mass [Kg], joint size [m]
    vec3 pos,vel,F;     // ball state [m],[m/s],[N]
    };

const int balls=3;
_ball ball[balls];
//---------------------------------------------------------------------------
void pendulum_update(double dt)
    {
    const float k=0.001;        // air friction coefficient should be function of object shape
    int i,j;
    float v;
    vec3 dir,nor,F0,F1;
    _ball *b,*b0,*b1;
    vec3 g=vec3(0.0,+9.81,0.0); // local gravity field
    dt*=10.0;                   // time multiplier for simulation speed ...

    // clear forces
    for (b=ball,i=0;i<balls;i++,b++) b->F=vec3(0.0,0.0,0.0);
    // compute driving force (F = m*acc) as gravity - air friction
    for (b=ball,i=0;i<balls;i++,b++) if (b->ix!=-2) b->F=(b->m*g) - b->vel*(k*length(b->vel));
    // apply bonds on F (in reverse order to avoid the need for recursion)
    for (b0=ball+balls-1,i=0;i<balls;i++,b0--)
     if (b0->ix>=0)
        {
        b1=ball+b0->ix;
        dir=normalize(b1->pos-b0->pos);     // unit direction vector of bond
        // remove radial part of acc
        F0=dir*dot(dir,b0->F);              // radial part of b0->F
        F1=dir*dot(dir,b1->F);              // radial part of b1->F
        F1=vec3(0.0,0.0,0.0);
        if (b0->ix!=-2){ b0->F+=F1-F0; }
        if (b1->ix!=-2){ b1->F+=F0-F1; }
        }

    // Newton/d'Lambert simulation vel = Integral((F/m)*dt)
    for (b=ball,i=0;i<balls;i++,b++) if (b->ix!=-2) b->vel+=(b->F/b->m)*dt;
    // apply bonds on vel
    for (b0=ball+balls-1,i=0;i<balls;i++,b0--)
     if (b0->ix>=0)
        {
        b1=ball+b0->ix;
        dir=normalize(b1->pos-b0->pos);     // unit direction vector of bond
        // transfer all vel to axial direction
        if (b0->ix!=-2)
            {
            v=length(b0->vel);
            if (v>1e-6)
                {
                nor=cross(b0->vel,dir);
                nor=cross(nor,dir);
                if (length(nor)>1e-6) b0->vel=-v*normalize(nor);
                }
            }
        }

    // Newton/d'Lambert simulation pos = Integral(vel*dt)
    for (b=ball,i=0;i<balls;i++,b++) if (b->ix!=2) b->pos+=b->vel*dt;
    // apply bonds on pos
    for (b0=ball+balls-1,i=0;i<balls;i++,b0--)
     if (b0->ix>=0)
        {
        b1=ball+b0->ix;
        dir=normalize(b1->pos-b0->pos);     // unit direction vector of bond
        // repair b0 distance to b1
        b0->pos=b1->pos-b0->l*dir;
        }
    }
//---------------------------------------------------------------------------
void pendulum_init(float xs,float ys)               // init pendulm to screen with xs,ys dimensions
    {
    Randomize();
    vec3 p,v=vec3(0.0,0.0,0.0);
    _ball *b=ball;
    // init balls
    p=vec3(0.5*xs,0.25*ys,0.0);
    b->ix=-2; b->m=0.5; b->r=0.015*xs; b->pos=p; b->vel=v; if (b->ix>=0) b->l=length(b->pos-ball[b->ix].pos); else b->l=0.0; b++; p.x+=0.15*xs;
    b->ix= 0; b->m=0.5; b->r=0.020*xs; b->pos=p; b->vel=v; if (b->ix>=0) b->l=length(b->pos-ball[b->ix].pos); else b->l=0.0; b++; p.x+=0.15*xs;
    b->ix= 1; b->m=0.5; b->r=0.020*xs; b->pos=p; b->vel=v; if (b->ix>=0) b->l=length(b->pos-ball[b->ix].pos); else b->l=0.0; b++;
    }
//---------------------------------------------------------------------------

I added mass m and change acceleration acc to force F parameters on per ball manner. I also added bond as index to which ball current ball is connected...

Here preview:

So I just precompute all the F ... then transfer radial parts of F along bonds to joined ball ... and only then apply the integration.