Wave equation with explicit Euler time integration


I am trying to solve the wave equation using explicit Euler time stepping (for the sake of simplicity) by substituting \dot{p}= q.
Thus, I am solving the set of equations summarized in the attached PDF-file.

For 1D, I am running the attached (executable) code, where at the Dirichlet boundary, the pressure p(t)=\sin/(2 \pi f t) with $f=1$Hz is prescribed.

Plots of the solution of the FEM simulation and of the extracted time signals at the Dirichlet BC at x=0 for p (which is prescribed) and q which is obtained from the simulation are attached.
The solution of the wave propagation, i.e. pressure field p is as expected. However, the FEM solution of the auxiliary quantity q is not. But when I am computing the difference quotient from p, I get the correct (target) solution for q.
Since the solution (of the primary quantity p) is right, I am wondering if I am using some of the NGSolve utilities wrong or if there is any other reason for that.

I also attached results of alternative excitation functions for the Dirichlet BC.
For the 2D case and also for an initial condition (instead of excitation by Dirichlet BC), I get the same behavior.