Problem with the numerical integration of the wave equation for an stretched string

Question

I am integrating the wave equation for an stretched string with an explicit finite difference method and I am getting an weird oscillation in the total energy (mechanical energy) of the string, whereas in theory it should be a constant of time. Here's a figure with the initial and final profiles (after 10.5 periods):

And here is a time series of the mechanical energy of the string

The code in python is here

# -*- coding: utf-8 -*-
#!/usr/bin/python
'''
waveEqFinDiff.py
Finite differences integration of 1d wave eq. of motion for a stretched string
'''
import numpy as np
import matplotlib.pyplot as plt

L = 1.0 # Length of the support
N = 256 # Number of grid points
dx = 1.0/(N-1) # smallest grid space
X = np.linspace(0, 1.0, N, endpoint = True)
k = 4*np.pi/L # wave constant
# period, time-step, time interval of integration
T = 2*np.pi/k

dt = T/1024
tMax = 10.5*T
Delta_t = np.arange (0, tMax, dt)
print "dx = %g, dt = %g, dt/dx = %g" % (dx, dt, dt/dx)
###########################################################################

def get_totEnergy(X, Y, dYdt, dt, dx):
    # get potential energy
    Y_x = (Y[2:] -Y[:-2])/(2*dx) #central difference 1st-order derivative
    Y_x = np.concatenate([[0.0], Y_x, [0.0]])
    # forward difference derivative
    # left border right-sided derivative approximation to second-order
    Y_x[0]= (4*Y[1]-Y[2])/(2*dx)
    # backward difference derivative
    # right border left-sided derivative approximation to second-order
    Y_x[-1]= (-4*Y[-2]+Y[-3])/(2*dx)
    E_pot= 0.5*np.sum(Y_x**2)*dx
    # get kinetic energy
    E_kin = 0.5*np.sum(dYdt**2)*dx
    E = E_kin+E_pot # mechanical energy
    return E
###########################################################################
'''
Integration of the wave equation for a stretched string using explicit 
finite differences method 
'''
def integWaveEq(Y, Delta_t, dx):
    Y_dt = Y
    totE = np.zeros(len(Delta_t))
    for idx, t in enumerate(Delta_t):
        #Y_xx
        Y_xx = (Y[2:]+Y[:-2]-2*Y[1:-1])/(dx**2)
        # boundary conditions:  Y[0, t]=Y[L, t]=0
        Y_xx = np.concatenate([[0.0], Y_xx, [0.0]])
        temp = Y
        # Time increment
        Y = 2*Y-Y_dt+dt*dt*Y_xx
        dYdt= (Y-Y_dt)/(2*dt)
        Y_dt = temp
        totE[idx] = get_totEnergy(X, Y, dYdt, dt, dx)
    return [Y, totE]
'''

'''
# Initial conditions of the string
Y = 0.1*np.sin(k*X)
plt.figure()
plt.plot(X, Y, 'r-', label = u"initial profile")
# Integrate the pde (wave equation)
Y, E_tot= integWaveEq(Y, Delta_t, dx)
# Check Von Neumann stability criterion
alpha = dt/dx
beta = 1.0-2*(alpha*np.sin(k*dx/2))**2
print "If |beta|<1, integration is stable"
print "beta=", beta
#plt.axes().set_aspect('equal')
plt.plot(X, Y, 'b-', label = u"final profile")
plt.legend()
plt.savefig("profiles.png")
plt.grid()
plt.figure()
plt.title("Total energy")
plt.plot(Delta_t, E_tot, 'b-')
plt.xlabel("t")
plt.ylabel("$E_{tot}$")
plt.grid()
plt.savefig("totEnergy.png")
plt.show()

I wonder how to improve the numerical integration, e.g. that could be verified by a reduction of the variation of the total energy, that should be conserved in time.