# Using NIntegrate and DiscretePlot to visualize pseudodifferential operators

In harmonic analysis, pseudodifferential operators are a way to generalize the notions of derivatives, through the use of Fourier transforms. The basic idea being,

1. Let $u(x)\in\mathcal{S}(\mathbb{R}^n)$ be a Schwartz function.
2. Take the Fourier transform to get $\hat{u}(\xi)\in\mathcal{S}(\mathbb{R}^n)$, since the Fourier transform is an isometry of $\mathcal{S}$ onto itself.
3. Multiply $\hat{u}(\xi)$ by a "symbol". That is, a function $P(x,\xi)\in C^{\infty}(\mathbb{R}^n)\times C^{\infty}(\mathbb{R}^n)$, with some complicated restrictions on how fast it and its derivatives can grow.
4. Take the inverse Fourier transform to get back a function of $x$, usually denoted as $(P(x,D)u)(x)$.

I know Mathematica has the F.T. built in already, but I'd like to be able to use NIntegrate methods for cases when the F.T. operation fails. For example, the first example I'm trying to work on involves $u(x)$ as a compactly supported function, u[x_]:=Exp[1]*Exp[-1/(1 - x^2)] for $|x|<1$ and $u\equiv0$ for $|x|\geq 0$ (defined piecewise in Mathematica.)

I'm also not (yet) looking into examples outside of $\mathbb{R}^1$.

When written out, as two integrals in $\mathbb{R}^1$, finding $(P(x,D)u)(x)$ is the same as taking $$(P(x,D)u)(x)=(2\pi)^{-1}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} e^{i(x-y)\xi}P(x,\xi)u(y)d\xi dy.$$ Notice that at least in my first example, the bounds on the $dy$ integral can be changed to just $\int_{-1}^{1}\cdots dy$ because of the compact support on $u(x)$.

I define a symbol function, p[x_,ξ_], and naively defined the function Pu[x_],

Pu[x_]:=(1/(2Pi))*NIntegrate[p[x,ξ]*NIntegrate[Exp[I*(x-y)*ξ]*u[y],{y,-1,1}],{ξ,-10,10}],

as two nested NIntegrate functions. It was my intention to use DiscretePlot to sample a number of values of $x$, and get a rough idea of what Pu[x] looks like (plotting the real and imaginary parts separately).

This works..."okay"... Right away it spits out these errors:

"The integrand ... has evaluated to non-numerical values for all sampling points in the region with boundaries {{-1,1}}"

"General::stop: Further output of NIntegrate::will be suppressed during this calculation"

"NIntegrate::ncvb : NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in recursive bisections in ξ {ξ}={4.14}. NIntegrate obtained -1.05471*10^-15+1.13798*10^-10 I and 2.32456*10^-11 for the integral and error estimates."

And after a considerable amount of time, it finally prints out what I'm assuming is the correct graph. But, considering how much time it takes to get even one plotting, I cannot really play around with it at all. Changing the plot range; integrating $\xi$ over a larger interval to get a better approximation; trying new symbols, $P(x,\xi)$; etc... these are all things I would love to tweak, but simply am unable to due to the computation times.

Does anyone here have any suggestions on how I could streamline this process? I feel like this is a great opportunity to visualize results from some currently relevant mathematics. I haven't seen any pictures out there, but I'm having trouble contributing with my horribly inefficient code.

Also, naturally, if we can get the process working in $\mathbb{R}^1$, then the next step would be to get it working in $\mathbb{R}^2$.

edit: Thanks to the suggestions, I now have pared down my function into its real and imaginary parts, each defined with only one NIntegrate:

realPu[x_]:=(1/(2Pi))*NIntegrate[p[x,ξ]*Cos[(x-y)*ξ]*u[y],{y,-1,1},{ξ,-30,30}] imgPu[x_]:=(1/(2Pi))*NIntegrate[p[x,ξ]*Sin[(x-y)*ξ]*u[y],{y,-1,1},{ξ,-30,30}]

Overall though, I don't see any noticeable improvement. Maybe someone has a completely new suggestion? I wonder if there is a totally different approach I should be taking instead.

edit 2: So I chose my symbol to be $P(x,\xi)=\log(1+\xi^2)$ (which corresponds to the operator $\log(1-\Delta)$) and let it run it's course. Here is the graph of the original bump function and its pseudoderivative:

edit 3: And here is the same $\Psi DO$ applied to a small perturbation of the original bump function:

• How is P[x_,ξ_] defined? (Give an example, I mean.) Sep 4, 2015 at 5:46
• Well, there could be many different symbols. The first one I tried was $P(x,\xi)=\xi$, as this should have been equivalent to simply taking one derivative classically. The limited data I was able to see seemed to roughly agree with this. The next symbol I'd like to try would be $P(x,\xi)=\log(1+\xi^2)$, or $P(x,\xi)=\cos(x)\log(1+\xi^2)$. Sep 4, 2015 at 5:49
• A couple of comments: You don't need to nest the two NIntegrates: the function supports multi-dimensional integration. Second: could be a copy-and-paste error, but imaginary unit is capital I in Mathematica. Third: it is good MMA practice to not capitalize user-defined functions so as not to interfered with built-in MMA symbols and functions (which all start with a capital letter). Sep 4, 2015 at 5:52
• I am not sure why, but I could speed up integration by splitting the Exp in Sin and Cos and writing it in two integrals, one real one imaginary. I guess it depends on how NIntegrate handles oscillatory integrands, while intuitively I would not have expected this to change anything. Sep 4, 2015 at 7:48
• Thanks. I now have two functions, each of which is defined using only one NIntegrate: realPu[x_] and imgPu[x_]. Running them with the Sin and Cos separated doesn't seem particularly faster though. Sep 4, 2015 at 21:06