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 quantitiy 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.

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