# Implementing velocity verlet for harmonic oscillator in C gives error wrt conservation of energy

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;
t=0;
ft = 10.;
dt = .01;
x0 = 0;
v0 = 1;
m = 1;
k= 1;
omega = .1;

for(t=0;t<=ft;t+=dt)
{

if(t==0)
{
x = position(x0, omega, v0,t);
v = velocity (x0, omega, v0, t);
a = acceleration(x0, omega, v0, t);
}

else
{
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;
}

• How can k=1 and m=1, but omega=.1? Doesn't $\omega=\sqrt{\frac{k}{m}}$ – Tyberius Apr 19 '17 at 17:03
• @Tyberius hey, you're correct, but for this case those values are arbitrary – L L Apr 19 '17 at 17:39
• 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. – Tyberius Apr 19 '17 at 17:41

## 1 Answer

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.

• Thank you so much!! That was exactly it. Hope you have a great day! – L L Jun 22 '16 at 18:56