Assume that $p(x)$ and $f(x)$ are sufficiently smooth functions and $D\equiv \frac{d}{dx}$. My question is concerned with the discretization of $p(x+D)f(x)$.
As an example, let $p(x)=x^{2}+2x$. Then we have $p(x+D)f(x)=(D^{2}+4xD+x^{2}+2x)f(x)$ and this can be discretized using the finite difference method.
Question: what can be done if $p(x)$ is not a polynomial (e.g., $p(x)=sin^{2}(x)$) or it is only known numerically?
My immediate thought is to approximate $p(x)$ with a polynomial and then use the technique described above. I however do not find it very elegant, and most importantly the error in this approach will be affected by the accuracy of the polynomial approximation. I appreciate it if you could suggest an alternative approach. My hunch is that maybe Fourier/Laplace transform is the way to go but do not know how.
Update: I have found an interesting property of $p(x+\frac{d}{dx})$ operator:
Let $T_{x}=(x+\frac{d}{dx})$, $T_{\omega}=i(\omega+\frac{d}{d\omega})$ and $F\{f\}=\hat f(\omega)$, using Fourier transform we have:
$$F\{T_{x}f\}=F\{xf+\frac{df}{dx}\}=F\{xf\}+F\{\frac{df}{dx}\}=i\frac{d}{d\omega} \hat f(\omega)+i\omega\hat f(w)=i(\omega+\frac{d}{d\omega}) \hat f(\omega)=T_{\omega} \hat f(\omega)$$.
Additionally, let us assume $g(x)=T_{x}f$, then $F\{T_{x}^{2}f\}=F\{T_{x} T_{x} f\}=F\{T_{x}g(x)\}=T_{\omega} \hat g(\omega)=T_{\omega} T_{\omega} \hat f(\omega)=T_{\omega}^{2}\hat f(\omega)$. If we repeat the above process, we have $F\{T_{x}^{k}f\}=T_{\omega}^{k} \hat f(\omega)$, therefore, if we assume $p(x)=\sum_{k=0}^\infty p_{k}x^k$, we have:
$$F\{p(x+\frac{d}{x})f(x)\}=F\{\sum_{k=0}^\infty p_{k}(x+\frac{d}{dx})^k f(x)\}=F\{\sum_{k=0}^\infty p_{k}T_{x}^k f(x)\}=\sum_{k=0}^\infty p_{k} T_{\omega}^k \hat f(\omega)= p(T_{\omega}) \hat f(\omega)= p(i(\omega+\frac{d}{d\omega})) \hat f(\omega).$$
The operator $p$ remained unchanged under Fourier transform (except for the multiplication with $i$).
