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.
We usually define $p=\hbar k$. In this post, I'll take $\hbar =1$ (see natural units) to simplify the notation. This means that $k=p$. This is not necessary, but IMHO this makes the arguments more transparent (besides the fact that natural units are indeed more natural once you get used to them). I usually recommend every learning physicist to try to use natural units whenever they can. If you don't like to, you can exchange $k\leftrightarrow p/\hbar$ in my answer.
In general, the expectation value of the momentum is given by
$$
\langle k\rangle= \int_{-\infty}^{+\infty} \mathrm dk\ k\ |\phi(k)|^2 \tag{1}
$$
which means that $|\phi(k)|^2$ has to decrease faster than $1/k^3$ as $k\to\infty$ for $\langle k\rangle$ to be well defined. In fact, we also want $\langle k^2\rangle$ to be well defined, which means we usually want $\phi(k)$ to be $\mathcal O(k^{-2})$ as $k\to\infty$.
If $\phi(k)$ has a large spread (it doesn't decrease fast enough), then the momentum is not well-defined, because the integrals that define momentum are divergent.
Why Does $\phi(k)$ have such a strong influence on momentum?
Well, $\phi(k)$ is the wave-function in momentum space, so all the information about the momentum of the system is contained in $\phi(k)$.
And how does a small spread in $\phi(k)$ correspond to a more well defined position?
Note that $\phi(k)$ is the Fourier Transform of $\Psi$, so that this statement is actually equivalent to the uncertainty principle of the Fourier Transform.
UPDATE
Here we prove that both expressions $(1)$ and $(2)$ are equivalent:
$$
\langle k\rangle=-i\int_{-\infty}^{+\infty} \mathrm dx \ \Psi^*\partial_x\Psi \tag{2}
$$
The proof is non-trivial, but in the end is just the differentiation property of the Fourier Transform. To prove the equivalence, we'll make use of the known fact
$$
\int_{-\infty}^{+\infty} \mathrm dx\ \mathrm e^{iqx}=2\pi\delta(q) \tag{3}
$$
for any $q\in\mathbb R$. Here, $\delta(q)$ is the Dirac delta function.
First, note that
$$
\Psi=\frac{1}{\sqrt{2\pi}}\int \mathrm dk\ \phi(k)\;\mathrm e^{i(kx-\omega t)} \tag{4}
$$
so that
$$
\partial_x\Psi=\frac{i}{\sqrt{2\pi}}\int \mathrm dk\ k\;\phi(k)\;\mathrm e^{i(kx-\omega t)} \tag{5}
$$
If we make the change of variable $k\to k_1$ in $(4)$ and $k\to k_2$ in $(5)$, and plug these into $(2)$ we get
$$
\begin{align}
\langle k\rangle&=\frac{1}{2\pi}\int \mathrm dx \int\mathrm dk_1\int\mathrm dk_2\ \overbrace{\phi^*(k_1)\;\mathrm e^{-i(k_1x-\omega t)}}^{\Psi^*}\ \overbrace{k_2\;\phi(k_2)\;\mathrm e^{i(k_2x-\omega t)}}^{\partial_x\Psi}=\\
&=\frac{1}{2\pi}\int\mathrm dk_1\int\mathrm dk_2\ k_2\ \phi^*(k_1)\;\phi(k_2)\int \mathrm dx\ \mathrm e^{i(k_2-k_1)x}
\end{align}
$$
Next, use $(3)$ to simplify the $x$ integral, where $q\equiv k_2-k_1$:
$$
\begin{align}
\langle k\rangle&=\int\mathrm dk_1\int\mathrm dk_2\ k_2\ \phi^*(k_1)\;\phi(k_2)\;\delta(k_2-k_1)=\\
&=\int\mathrm dk_1 k_1\ |\phi(k_1)|^2
\end{align}
$$
which is just $(1)$ upon the change of variables $k_1\to k$. As you can see, both expressions are equivalent, which means that you can use whichever you like the most.
Best Answer
These relations are found in every book on QM, but the usual notation is $$ X|x\rangle=x|x\rangle $$ and $$ P|p\rangle=p|p\rangle $$
To go from these equations to the ones you've written, you just have to project them into the position basis $|x'\rangle$ (and use $\langle x'|x\rangle=\delta(x-x')$ and $\langle x'|p\rangle\sim\exp[ipx]$).
Edit Important: see ACuriousMind's comment below.