  Linear wave equations sourced by a Dirac delta distribution $\delta(x)$ and
its derivative(s) can serve as a model for many different phenomena. We
describe a discontinuous Galerkin (DG) method to numerically solve such
equations with source terms proportional to $\partial^n \delta /\partial x^n$.
Despite the presence of singular source terms, which imply discontinuous or
potentially singular solutions, our DG method achieves global spectral accuracy
even at the source's location. Our DG method is developed for the wave equation
written in fully first-order form. The first-order reduction is carried out
using a distributional auxiliary variable that removes some of the source
term's singular behavior. While this is helpful numerically, it gives rise to a
distributional constraint. We show that a time-independent spurious solution
can develop if the initial constraint violation is proportional to $\delta(x)$.
Numerical experiments verify this behavior and our scheme's convergence
properties by comparing against exact solutions.
