A partial Fourier transform method for a class of hypoelliptic Kolmogorov equations

We consider hypoelliptic Kolmogorov equations in $n+1$ spatial dimensions, with $n\geq 1$, where the differential operator in the first $n$ spatial variables featuring in the equation is second-order elliptic, and with respect to the $(n+1)$st spatial variable the equation contains a pure transport term only and is therefore first-order hyperbolic. If the two differential operators, in the first $n$ and in the $(n+1)$st co-ordinate directions, do not commute, we benefit from hypoelliptic regularization in time, and the solution for $t>0$ is smooth even for a Dirac initial datum prescribed at $t=0$. We study specifically the case where the coefficients depend only on the first $n$ variables. In that case, a Fourier transform in the last variable and standard central finite difference approximation in the other variables can be applied for the numerical solution. We prove second-order convergence in the spatial mesh size for the model hypoelliptic equation $\frac{\partial u}{\partial t} + x \frac{\partial u}{\partial y} = \frac{\partial^2 u}{\partial x^2}$ subject to the initial condition $u(x,y,0) = \delta (x) \delta (y)$, with $(x,y) \in \mathbb{R} \times\mathbb{R}$ and $t>0$, proposed by Kolmogorov, and for an extension with $n=2$. We also demonstrate exponential convergence of an approximation of the inverse Fourier transform based on the trapezium rule. Lastly, we apply the method to a PDE arising in mathematical finance, which models the distribution of the hedging error under a mis-specified derivative pricing model

Near-optimal perfectly matched layers for indefinite Helmholtz problems

A new construction of an absorbing boundary condition for indefinite Helmholtz problems on unbounded domains is presented. This construction is based on a near-best uniform rational interpolant of the inverse square root function on the union of a negative and positive real interval, designed with the help of a classical result by Zolotarev. Using Krein's interpretation of a Stieltjes continued fraction, this interpolant can be converted into a three-term finite difference discretization of a perfectly matched layer (PML) which converges exponentially fast in the number of grid points. The convergence rate is asymptotically optimal for both propagative and evanescent wave modes. Several numerical experiments and illustrations are included.Comment: Accepted for publication in SIAM Review. To appear 201

Preconditioners for the spectral multigrid method

The systems of algebraic equations which arise from spectral discretizations of elliptic equations are full and direct solutions of them are rarely feasible. Iterative methods are an attractive alternative because Fourier transform techniques enable the discrete matrix-vector products to be computed with nearly the same efficiency as is possible for corresponding but sparse finite difference discretizations. For realistic Dirichlet problems preconditioning is essential for acceptable convergence rates. A brief description of Chebyshev spectral approximations and spectral multigrid methods for elliptic problems is given. A survey of preconditioners for Dirichlet problems based on second-order finite difference methods is made. New preconditioning techniques based on higher order finite differences and on the spectral matrix itself are presented. The preconditioners are analyzed in terms of their spectra and numerical examples are presented

Regular polynomial interpolation and approximation of global solutions of linear partial differential equations

We consider regular polynomial interpolation algorithms on recursively defined sets of interpolation points which approximate global solutions of arbitrary well-posed systems of linear partial differential equations. Convergence of the 'limit' of the recursively constructed family of polynomials to the solution and error estimates are obtained from a priori estimates for some standard classes of linear partial differential equations, i.e. elliptic and hyperbolic equations. Another variation of the algorithm allows to construct polynomial interpolations which preserve systems of linear partial differential equations at the interpolation points. We show how this can be applied in order to compute higher order terms of WKB-approximations of fundamental solutions of a large class of linear parabolic equations. The error estimates are sensitive to the regularity of the solution. Our method is compatible with recent developments for solution of higher dimensional partial differential equations, i.e. (adaptive) sparse grids, and weighted Monte-Carlo, and has obvious applications to mathematical finance and physics.Comment: 28 page

A fast immersed boundary method for external incompressible viscous flows using lattice Green's functions

A new parallel, computationally efficient immersed boundary method for solving three-dimensional, viscous, incompressible flows on unbounded domains is presented. Immersed surfaces with prescribed motions are generated using the interpolation and regularization operators obtained from the discrete delta function approach of the original (Peskin's) immersed boundary method. Unlike Peskin's method, boundary forces are regarded as Lagrange multipliers that are used to satisfy the no-slip condition. The incompressible Navier-Stokes equations are discretized on an unbounded staggered Cartesian grid and are solved in a finite number of operations using lattice Green's function techniques. These techniques are used to automatically enforce the natural free-space boundary conditions and to implement a novel block-wise adaptive grid that significantly reduces the run-time cost of solutions by limiting operations to grid cells in the immediate vicinity and near-wake region of the immersed surface. These techniques also enable the construction of practical discrete viscous integrating factors that are used in combination with specialized half-explicit Runge-Kutta schemes to accurately and efficiently solve the differential algebraic equations describing the discrete momentum equation, incompressibility constraint, and no-slip constraint. Linear systems of equations resulting from the time integration scheme are efficiently solved using an approximation-free nested projection technique. The algebraic properties of the discrete operators are used to reduce projection steps to simple discrete elliptic problems, e.g. discrete Poisson problems, that are compatible with recent parallel fast multipole methods for difference equations. Numerical experiments on low-aspect-ratio flat plates and spheres at Reynolds numbers up to 3,700 are used to verify the accuracy and physical fidelity of the formulation.Comment: 32 pages, 9 figures; preprint submitted to Journal of Computational Physic