I have a Wiener process in 3 dimensions, with time-independent drift and diffusion, ie x \in \mathbb{R}^3,
dx_t = \mu(x) dt + \sigma(x) dw
where \sigma(x) is diagonal. I also have an initial distribution/density p_0(x) at t = 0. Importantly, the domain is infinite (ie x is not bounded).
I would like to obtain an approximation of the distributions, eg densities p_t(0). I have already simulated this and everything is nice and smooth, so I am wondering if I could be using spectral methods for this.
Reading up about the problem, almost all texts that I have seen recommend pseudospectral methods where they convert to an ODE and solve that. But I would like to understand why one cannot treat the whole domain (ie \mathbb{R}^3 \times [0, T]) with spectral methods. The only package I have seen that seems to do this is ApproxFun.jl.
Is there a numerical problem with using a spectral representation for every dimension, including time, that the approach of Olver and Townsend somehow avoids? Or is the approximation is OK, I could be using eg a rationally transformed Chebyshev basis as well?
(Note: answers can assume that I am clueless about numerical solutions to PDEs and recommend nice intro texts. I have Boyd (2000)).
If you use tensor product Chebyshev basis then 3d already has very big matrices. So 4d will be too big to be practical. (And I don’t think it actually works in ApproxFun.jl. )
That said, if your problem is on ℝ^3 and if your variable coefficients are rotationally invariant then you could use Zernike polynomials in a ball tensored with Chebyshev. This would decouple across spherical harmonic modes and hence your actual systems will be 2D instead of 4D. Unfortunately we don’t have the code to do balls yet (just disks in MultivariateOrthogonalPolynomials.jl)
In principle, no problem. You can do the same thing in time as spatial dimension. PINNs do it this way.
But as mentioned, if there is a way to do time stepping, that will save your resources since you don’t need to store all the information in time axis.
Thanks for all the answers. I have experimented with various approaches and spectral approximations are definitely viable for my purposes.
The only difficulty I am running into is ensuring that the approximated density is positive. Eg if I approximate
p(x, t) = \sum_j c_j \cdot \varphi(x, t)
where p(x, t) \ge 0 and \int p(x, t) dt = 1, ie a proper density function that satisfies the relevant PDE, the linear form makes enforcing the integral really easy, but does not rule out p(x, t) < 0.
At this point I do not know if this is a problem for me in practice, as generally when the density is below 0 it is really tiny in magnitude, so I could just treat it as numerical error, or use \max(0, p(x, t)) instead. But if anyone has suggestions on how this is usually handled, I would appreciate them.