Improved Euler Method in Simple Harmonic Oscillator

Question

I have written a C code using the improved Euler method to determine the position, velocity and energy of the oscillator at regular time intervals. However, I run into a problem that the energy of the oscillator is decreasing, though there are no dissipation terms. I think this is particularly related with the way I update my position and velocity variables and would like to get your help on the matter. My code is as follows:

//Compilation and run
//gcc oscillatorimprovedEuler.c -lm -o oscillatorimprovedEuler && ./oscillatorimprovedEuler
#include <stdio.h> 
#include <math.h>  

// The global constans are defined in the following way (having the constant value througout the program
#define m 1.0                 // kg
#define k 1.0                 // kg/sec^2
#define h 0.1                // sec This is the time step
#define N 201 // Number of time steps

int main(void)
{
    // We avoid using arrays this time
    double x = 0, xint = 0;
    double v = 5, vint = 0; // Just like the previous case
    double t = 0;
    double E = (m * v * v + k * x * x) / 2.0; // This is the energy in units of Joules

    FILE *fp = fopen("oscillatorimprovedEuler.dat", "w+");
    int i = 0;
    for(i = 0; i < N ; i++)
    {
        fprintf(fp, "%f \t %f \t %f \t %f \n", x, v, E, t);
        xint = x + (h) * v;
        vint = v - (h) * k * x / m;
        v = v - (h) * ((k * x / m) + (k * xint / m)) / 2.0;
        x = x + (h) * (v + vint) / 2.0;
        E = (m * v * v + k * x * x) / 2.0;
        t += h;
    }
    fclose(fp);
    return 0;
}

There may be a very slight point I miss so I would be grateful if you can point it out. I appreciate your help.

Answer 1

So I figured out with the aid of math.stackexchange that the problem was related with updating the position and velocity earlier than the time they should be updated and more intermediate variables were needed. The now working code is below:

//Compilation and run
//gcc oscillatorimprovedEuler.c -lm -o oscillatorimprovedEuler && ./oscillatorimprovedEuler
#include <stdio.h> 
#include <math.h>  

// The global constans are defined in the following way (having the constant value througout the program
#define m 1.0                // kg
#define k 1.0                // kg/sec^2
#define h 0.1                // sec This is the time step
#define N 200                // Number of time steps

int main(void)
{
    // We need to define this many variables to avoid early updating the position and velocity
    double x = 0.0, xpre = 0, xcor = 0;
    double v = 5.0, vpre = 0, vcor = 0; // Just like the previous case
    double t = 0;
    double E = (m * v * v + k * x * x) / 2.0; // This is the energy in units of Joules

    FILE *fp = fopen("oscillatorimprovedEuler.dat", "w+");
    int i = 0;
    for(i = 0; i < N ; i++)
    {
        if (i == 0)
        {
            fprintf(fp, "%f \t %f \t %f \t %f \n", x, v, E, t);
        }
        xpre = x + (h) * v;
        vpre = v - (h) * k * x / m;
        vcor = v - (h) * ((k * x / m) + (k * xpre / m)) / 2.0;
        xcor = x + (h) * (v + vpre) / 2.0;
        E = (m * vcor * vcor + k * xcor * xcor) / 2.0;
        t += h;
        fprintf(fp, "%f \t %f \t %f \t %f \n", xcor, vcor, E, t);
        x = xcor, v = vcor;
    }
    fclose(fp);
    return 0;
}