I'm attempting to implement velocity Verlet with a harmonic oscillator in C, but I have some errors:

  1. Energy is not being conserved
  2. Energy oscillates much more than I would expect

I believe the error has to do with my position function or definitions, as the PE is blowing up much faster than the KE is going down, but I have not been able to determine the error. It's also possible that the KE is not going down as fast as it should, I suppose. I guess it should be mentioned that I'm new to programming in C.

#include <stdio.H>
#include <math.H>

double position(double x0, double omega, double v0, double t)
    return x0 * cos( omega * t) + (v0/omega) * sin( omega * t);

double velocity(double x0, double omega, double v0, double t)
    return -x0 * omega * sin(omega * t) + v0 * cos( omega * t);

double acceleration(double x0, double omega, double v0, double t)
    return -x0 * omega*omega * cos( omega * t) - v0 * omega * sin( omega * t);

int main()
    double x0, v0, ao, x, xu, v, a, pe, ke, te, ft, t, m , dt, k, omega;
    ft = 10.;
    dt = .01;
    x0 = 0;
    v0 = 1;
    m = 1;
    k= 1;
    omega = .1;


            x = position(x0, omega, v0,t);
            v = velocity (x0, omega, v0, t);
            a = acceleration(x0, omega, v0, t);

            a = (-k * x)/m;
            v = v + (0.5*a*dt);
            x = x + (v * dt);
            ao = acceleration(x, omega, v, t);
            v = v + (0.5*ao*dt);

        ke= 0.5 * m * v * v;
        pe = 0.5 * k * x * x;
        te = pe + ke;
        //printf("position at time %f = %f\n",t,x);
        //printf("velocity at time %f = %f\n",t,v);
        //printf("acceleration at time %f = %f\n",t,a);*/
        printf("energy at time %f = K: %f; V: %f; Total: %f\n\n",t,ke,pe,te);


    return 0;
  • $\begingroup$ How can k=1 and m=1, but omega=.1? Doesn't $\omega=\sqrt{\frac{k}{m}}$ $\endgroup$ – Tyberius Apr 19 '17 at 17:03
  • $\begingroup$ @Tyberius hey, you're correct, but for this case those values are arbitrary $\endgroup$ – Luciano Apr 19 '17 at 17:39
  • $\begingroup$ Well sure you can choose two of them arbitrarily, but the third one is determined by your choice of the other two and, since you use all three of the values throughout the code, they need to consistent. $\endgroup$ – Tyberius Apr 19 '17 at 17:41

Your second acceleration ao in the else statement does not make sense. Since it assumes that the calculated position and velocity are the initial conditions at $t=0$ and then you evaluate it at the current time.

If you change the expression for ao to the same as you used for a, but using the updated position, then you should get the expected results.

| cite | improve this answer | |
  • $\begingroup$ Thank you so much!! That was exactly it. Hope you have a great day! $\endgroup$ – Luciano Jun 22 '16 at 18:56

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.