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 following set of equations:

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

The solution of the wave propagation, i.e. pressure field p is as expected (see following figures). 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.

Here the solution of the FEM simulation:

And the extracted time signals at the Dirichlet BC at x=0 for p (which is prescribed) and q which is obtained from the simulation:

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.