[Math] Improved Euler method for second order ODE

numerical methodsordinary differential equations

I am trying to solve the simple harmonic oscillator problem with various Euler methods. Having managed to solve it with simple and modified Euler methods now I am trying to solve it with the improved Euler method which is a bit like Runge-Kutta. My concern is that when I plot the graph of energy vs. time it oscillates though it also increases at the same time. Since x(position) is actually the integral of v(velocity) is it possible to solve it numerically?

The equations I have used are as follows:

$$ \frac{dv}{dt}= a= \frac{F}{m}=-\frac{kx}{m}$$

and the following first-order equation for the position:

$$ \frac{dx}{dt}= v $$

EDIT: The initial conditions and the constants are
$$x(0) = 0, v(0) = 5, \ \ k = m = 1,\ h=0.1, \ N=201 $$

Yet my concern is that these equations are not independent and I think the improved Euler method cannot be applied. Do you have any suggestions on how to approach the problem? I have also added the resulting graph to illustrate my point. I have updated position and velocity as follows and in the depicted order:
$$ x_{0}=x \\ x = x + hv \\ v = v – h(\frac{kx}{m} + \frac{kx_{0}}{m})/2 \\ E = (mv^2+kx^2)/2$$

The resulting graph after the data is plotted.

There must be a mistake in this implementation and I would be grateful if someone point it out or makes a suggestion.

Thanks in advance

EDIT: Latest C Implementation

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

Best Answer

It works for me (using this method):

solution

This is the spreadsheet I used (energy column not shown):

spreadsheet

The prediction is plain Euler: C2 = A2+B2*0.1, D2 = B2+0.1*(-A2).

The correction is trapezoidal: A3 = A2+0.1/2*(B2+D2), B3 = B2+0.1/2*(-A2-C2)

And so on. The way I typed the formulas in, the prediction on line 2 refers to the same time moment as row 3, which perhaps isn't logical. It could be placed one row below. But this is a matter of spreadsheet style. The point is, the method works.

The energy does grow slowly: at the end of the plot shown above, it is $12.6$ versus initial $12.5$. But this is perfectly normal for a numerical method.

Related Question