Unexpected noise in the solution of wave equation

Hello,

I'm solving the wave equations in 2D using the Newmark method for the time integration and something is wrong with my solution: it contains unexpected noise.

I'm solving it with respect to the acceleration and use the initial condition u0 on the displacement to find the initial condition a0 of the acceleration.

Can anyone spot where is the error?

from dolfin import *
from mshr import *
import numpy  as np




# Number of bumps of the waves
m_ = 1
n_ = 2
# Velocity
c = Constant(5.) 


# Order of the Finite Elements 
pdim = 1

# Create mesh
xlim = 2
ylim = 4
size = 10*3.1
domain = Rectangle(Point(0., 0.), Point(xlim, ylim))
mesh = generate_mesh(domain,size)

# Newmark Constants
beta_ = 0.25
gamma_ = 0.5

# Time Constants
dt =  0.001
T = 0.5
t = 0.


# Boundaries
class Horizontal(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and ((abs(x[1] - 0.) < DOLFIN_EPS) or (abs(x[1] - ylim) < DOLFIN_EPS))

class Vertical(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and ((abs(x[0] - 0.) < DOLFIN_EPS) or (abs(x[0] - xlim) < DOLFIN_EPS))


boundaries = FacetFunction("size_t", mesh, 0) 

Horizontal().mark(boundaries, 1) 
Vertical().mark(boundaries, 2) 


# Define Function space, trial and test functions
V = FunctionSpace(mesh, 'Lagrange', pdim)
a, phi = TrialFunction(V), TestFunction(V)

# Define Dirichlet BC
bcH = DirichletBC(V, Constant(0.), boundaries, 1)
bcV = DirichletBC(V, Constant(0.), boundaries, 2)
bcs = [bcH, bcV]


# Fields from previous time step (displacement, velocity, acceleration)
u0_expr = 'sin(m_*pi*x[0]/xlim)*sin(n_*pi*x[1]/ylim)'
u0IC = Expression(u0_expr, m_=m_, n_=n_, xlim=xlim, ylim=ylim, degree = pdim)
v0IC = Constant(0.)

u0 = interpolate(u0IC, V)
v0 = interpolate(v0IC, V)


# Blocks definition
def k_block(value):
    return inner(grad(value),grad(phi)) * dx 

m = pow(c,-2) * a * phi * dx 

# Assemble left hand side 
A = assemble(m + beta_*pow(dt,2)*k_block(a))
# Assemble right hand side
b0 = assemble(-k_block(u0))
# Impose BCs
[bc.apply(A) for bc in bcs]  

# Compute a0
a0 = Function(V)
solve (A,a0.vector(),b0)


# Save initial displacement (actually pressure here) and acceleration  
vtk_file_a = File("solution/p"+str(pdim)+"_"+str(size)+"_pressure_0acc.pvd")
vtk_file_a << a0