Commuting momentum and position operators? 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.)
 A: 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.   
A: 
Is there some requirement that observable operators should map C∞ functions to C∞ functions? In any case, I'd be grateful if someone could tell me why the above is ruled out.

By defining the position operator in the way you've indicated, you would introduce $artificial$ lattice of preferred values of coordinates in space $x_n = (n+1/2)\delta$ with (artificial) spacing $\delta$ and which would be stationary in arbitrary inertial frame. As you say, the new position and momentum operator violate canonical commutation relation, which is viewed as one of the basic features of operators in quantum theory.
By restricting the position operator to discrete eigenvalues $x_n$, one is restricting the set of possible positions to a discrete set. It is not clear why the other values of $x$ in between should be forbidden. The original Born  interpretation allows any value of $x$. Most naturally, possible configurations form a continuous set (imagine what would happen to your description if the frame of reference was shifted along $x$ by $\delta x <\delta$).
