Peter D. Drummond

Channa Hatharasinghe, Run Yan Teh, Jesse van Rhijn et al.

Spectral methods for partial differential equations (PDEs) with non-periodic boundary conditions arising in computational physics often use polynomial expansions on non-uniform grids. Here, we implement a Fourier method that employs fast trigonometric expansions on a uniform grid with non-periodic boundaries using fast discrete sine transforms (DST) or/and discrete cosine transforms (DCT) to solve parabolic PDEs. We implement this method in two ways: either using a Fourier spectral derivative or a Fourier interaction picture. Both methods can treat vector fields with a combination of Dirichlet and/or Neumann boundary conditions in one or more space dimensions. As examples, we use them to solve a variety of computational physics PDEs with analytical solutions, including the Peregrine solitary wave solution. For the 1D heat equation problem, our method with an interaction picture is accurate up to machine precision. Soluble examples of stochastic partial differential equation (SPDE) with non-periodic boundaries in one and two space dimensions, with physics and interdisciplinary applications are also treated. We compare the results obtained from these algorithms with publicly available solvers that use polynomial spectral methods, and study their relative performance and error scaling. Polynomial methods with non-uniform spatial grids have lower spatial discretization errors when the solutions change slowly in space, typically with large spatial grids. For problems with rapid spatial variation, Fourier methods can outperform polynomial expansions, owing to their smaller maximum space interval, and are generally faster due to the computational efficiency of discrete Fourier transform methods. We verified this by making a complexity analysis in which we studied the total error at the optimum combination of time and space steps for a given resource use.

Bogdan Opanchuk, Simon Kiesewetter, Peter D. Drummond

Stochastic equations play an important role in computational science, due to their ability to treat a wide variety of complex statistical problems. However, current algorithms are strongly limited by their sampling variance, which scales proportionate to 1/N_S for N_S samples. In this paper, we obtain a new class of variance reduction methods for treating stochastic equations, called parallel optimized sampling. The objective of parallel optimized sampling is to reduce the sampling variance in the observables of an ensemble of stochastic trajectories. This is achieved through calculating a finite set of observables - typically statistical moments - in parallel, and minimizing the errors compared to known values. The algorithm is both numerically efficient and unbiased. Importantly, it does not increase the errors in higher order moments, and generally reduces such errors as well. The same procedure is applied both to initial ensembles and to changes in a finite time-step. Results of these methods show that errors in initially optimized moments can be reduced to the machine precision level, typically around 10^(-16) in current hardware. For nonlinear stochastic equations, sampled moment errors during time-evolution are larger than this, due to error propagation effects. Even so, we provide evidence for error reductions of up to two orders of magnitude in a nonlinear equation example, for low order moments, which is a large practical benefit. The sampling variance typically scales as 1/N_S, but with the advantage of a very much smaller prefactor than for standard, non-optimized methods.