import math

def dx(v):
    return (v)

def dv(x, v, c, k, m):
    return (-c*v - k*x)/m

def e_mecanica(x, v, m , k):
    return (m*(v**2)+k*(x**2))/2

def Ralston(xi,vi,xin,vin,c,k,m,h):
    k1vt = dv(xin, vin, c, k, m)
    k1xt = dx(vin)
    k2vt = dv((xin + (3/4)*h), (vin + ((3/4)*k1vt*h)), c, k, m)
    k2xt = dx( (vin + (3/4)*k1vt*h))
    
    xin = xi + ((1/3)*k1xt + (2/3)*k2xt) * h
    vin = vi + ((1/3)*k1vt + (2/3) * k2vt) * h
    return(xin,vin)

m = 2
c = 0.5
k = 1
xi = 0
vi = 0.1
h = 0.1
t = 0
xin=xi
vin=vi

Em = e_mecanica(xi, vi, m , k) 
print ('Energia Mecanica Inicial= '+str(round(Em,4)))
Emnuevo = Em
t+=h
while (Emnuevo >= (Em/2)):
    xin,vin=Ralston(xi,vi,xin,vin,c,k,m,h)
    Emnuevo = e_mecanica(xin, vin, m , k) 
    print ('Energia Mecanica en '+ str(round(t,4))+' = '+str(round(Emnuevo,4)))
    xi = xin
    vi = vin
    t+=h
print('Tiempo: ' + str(round(t-h,2)))