I think this is a question about the mathematical axioms of quantum mechanics…
Consider the following operators $\tilde{x}$ and $\tilde{p}$ on Hilbert space $L^2(R)$, defined for fixed $\delta>0$:
$$
(\tilde{x}\psi)(x) = (n+1/2) \delta\ \psi(x)
$$
when $n\delta \leq x <(n+1)\delta$ for some integer $n$; and
$$
(\tilde{p}\psi)(x) = -i \hbar \frac{d\psi}{dx}(x)
$$
when $n\delta < x < (n+1) \delta$, and $(\tilde{p}\psi)(x) = 0$ when $x=n\delta$.
Both these operators are self-adjoint and they commute. Further, $\tilde{p}$ is identical to the usual momentum operator $\hat{p} = -i\hbar \frac{d}{dx}$ on the subspace of $L^2(R)$ on which $\hat{p}$ is defined, and thus has identical spectrum. Finally, by taking $\delta$ small we can make $\tilde{x}$ as close as we like to the usual position operator, $\hat{x}$. Thus naively it looks like $\tilde{x}$ and $\tilde{p}$ represent compatible observables which can be made arbitrarily close to $\hat{x}$ and $\hat{p}$.
I'm pretty sure that I've cheated somewhere, but I can't find anything in the axioms of quantum mechanics that rules out defining observables $\tilde{x}$ and $\tilde{p}$ as above (or at least, the axioms as they are usually presented in under/graduate physics texts). Is there some requirement that observable operators should map $C^\infty$ functions to $C^\infty$ functions? In any case, I'd be grateful if someone could tell me why the above is ruled out.
(As background, I am teaching an undergraduate class on quantum mechanics next semester and I wanted to explain to students that while it's not possible to make exact simultaneous measurement of position and momentum; it is possible to make an approximate simultaneous measurement. In other words, I wanted to find operators which approximate $\hat{x}$ and $\hat{p}$ on 'macroscopic' scales and which commute. This was before I looked in the literature and realised how technical constructing such approximate operators is.)
Best Answer
Concerning $\tilde{x}$, it is self-adjoint (the proof is easy) if it is defined on its natural domain $D(\tilde x)$. For $\gamma_n(x):= (n+1/2)\delta$ if $n \delta \leq x <(n+1)\delta$ and $n \in \mathbb Z$
$$D(\tilde{x}) = \left\{\psi\in L^2(\mathbb R) \:\left|\:\int_{\mathbb R} |\gamma_n(x)\psi(x)|^2 dx < +\infty \right.\right\}$$
The spectrum is $\sigma(\tilde{x})= \{\delta (n+1/2)\:|\: n \in \mathbb Z\}$ as expected.
The fact that your claimed approximated momentum operator $\tilde{p}$ makes any sense as it stands is questionable. The point is that $\tilde{p}$ defined this way is not self-adjoint but only symmetric, while observables must be self-adjoint in order to exploit the spectral calculus technology.
To obtain a symmetric operator (i.e. Hermitean and densely defined) I think that a possibility is to define its domaine $D(\tilde{p})$ as a space of $C^1$ functions (that vanish in $x_n= \delta n$ for all $n \in \mathbb Z$ with their first derivative for a reason I discuss shortly) also requiring that these functions rapidly vanish with the first derivative for $|x|\to \infty$. It should be possible to relax the smoothness condition using weak derivatives, but the situation does not seem to change remarkably.
Consider the anti-linear operator $(C\psi)(x):= \overline{\psi(-x)}$, the bar denoting the complex conjugation. It is norm preserving and $CC=I$, so it is a conjugation. It seems to me that $C\tilde{p}= \tilde{p}C$ with the domain I introduced above. Therefore, in view of a theorem due to von Neumann, $\tilde{p}$ admits some self-adjoint extension. A careful analysis of the defect indices of $\tilde{p}$ would classify these extensions. Therefore we are not authorized to say that $\tilde{p}$ has the meaning of an observable, we have to choose a self-adjoint extension of it. Unfortunately, the spectrum depends on this choice. It is by no means obvious what the spectrum of a self-adjoint extension of $\tilde{p}$ is.
What I am saying is that since we do not know the spectrum of the claimed approximated momentum we cannot say how (if) this claimed approximation works as soon as $\delta \to 0$.
It is possible to prove that if the domain is chosen as $C_0^\infty(\mathbb R)$ or ${\cal S}(\mathbb R)$ then there is a unique self-adjoint extension and coincides with the standard one. But we cannot make this choice in view of the form of $\tilde{x}$.
In fact, when computing $[\tilde{x}, \tilde{p}]=0$ we are assuming that the domain of $\tilde{p}$ is invariant under the action of $\tilde{x}$. This indeed happens if $D(\tilde{p})$ is made of the $C^1$ functions rapidly vanishing as $|x|\to \infty$ also vanishing at each $\delta n$ with their first derivative. The last condition eliminates the discontinuities introduced by the action of $\tilde{x}$ giving rise to a function in the same domain.
We conclude that $[\tilde{x}, \tilde{p}]=0$ holds on a domain which does not fix a unique self-adjoint extension of $\tilde{p}$ (and its spectrum consequently) as it would be if this domain were $C_0^\infty(\mathbb R)$ or ${\cal S}(\mathbb R)$.
An alternative definition of $\tilde{p}$ is obtained re-defining its domain $D(\tilde{p})$ by including the (vanishing sufficiently fast at infinity) functions $C^1$ in each open interval $(n\delta, (n+1)\delta)$ that are separately periodic in each such interval. The values these functions attain in the set of points $\delta n$ is irrelevant as this set has zero measure and $L^2$ does not care of zero measure set: We can safely assume that $\tilde{p}\psi$ vanishes thereon by definition as already assumed by the OP. This only concerns mathematics while, physically speaking, this assumption has devastating implications.
With this definition $\tilde{p}$ become essentially self-adjoint, i.e. it has a unique self adjoint extension. The spectrum should be $\sigma(\tilde{p}_{ext}) = \{ \frac{2\pi m}{\delta} \hbar\:|\: m \in \mathbb Z\}$, since in each said interval we find the standard momentum operator on a segment with periodic boundary conditions (to be completely sure I should check some conditions but I am reasonably confident)
With this definition you find $[\tilde{x}, \tilde{p}]=0$ on $D(\tilde{p})$.
The problem with this construction is that $\tilde{p}$ does not generate space displacements of the wavefunctions, $(e^{-ix_0\tilde{p}}\psi)(x)= \psi(x- x_0)$. And it does not seem that the situation improves as soon as $\delta \to 0$. Instead $\tilde{p}$ generates "periodic" displacements in each interval $(n\delta, (n+1)\delta)$ (as it were a circle). For $\delta \to 0$ these periodic displacements become denser and denser but they have nothing to do with geometric translations along $x$, as the ones generated by the true momentum operator.