Discontinuous collocation methods and gravitational self-force applications

Numerical simulations of extereme mass ratio inspirals, the mostimportant sources for the LISA detector, face several computational challenges. We present a new approach to evolving partial differential equations occurring in black hole perturbation theory and calculations of the self-force acting on point particles orbiting supermassive black holes. Such equations are distributionally sourced, and standard numerical methods, such as finite-difference or spectral methods, face difficulties associated with approximating discontinuous functions. However, in the self-force problem we typically have access to full a-priori information about the local structure of the discontinuity at the particle. Using this information, we show that high-order accuracy can be recovered by adding to the Lagrange interpolation formula a linear combination of certain jump amplitudes. We construct discontinuous spatial and temporal discretizations by operating on the corrected Lagrange formula. In a method-of-lines framework, this provides a simple and efficient method of solving time-dependent partial differential equations, without loss of accuracy near moving singularities or discontinuities. This method is well-suited for the problem of time-domain reconstruction of the metric perturbation via the Teukolsky or Regge-Wheeler-Zerilli formalisms. Parallel implementations on modern CPU and GPU architectures are discussed.Comment: 29 pages, 5 figure

Studies in numerical quadrature

Various types of quadrature formulae for oscillatory integrals are studied with a view to improving the accuracy of existing techniques. Concentration is directed towards the production of practical algorithms which facilitate the efficient evaluation of integrals of this type arising in applications. [Continues.

A formula for optimal integration in H2

AbstractThe weights âj of the optimal integration formula Q̂ = Σjâjf(zj) in H2 for given integration points zj are the exact integrals of the cardinal functions in the corresponding formula for optimal evaluation. By writing these cardinal functions as sums of their principal values, we very easily obtain a closed formula for the weights. In the case of real zj's, this formula makes explicit a series formula of Wilf. We compare numerically the accuracy of the optimal formula with that of some well-known integration formulae. For points equidistant on a circle of radius r, the formula allows an alternate derivation of a formula obtained by Golomb. We give also the barycentric formula for optimal evaluation with these points, as well as an experimentally stable sequence of radii r for integrating with an increasing number of points

Asymptotic expansions and fast computation of oscillatory Hilbert transforms

In this paper, we study the asymptotics and fast computation of the one-sided oscillatory Hilbert transforms of the form $$H^{+}(f(t)e^{i\omega t})(x)=-int_{0}^{\infty}e^{i\omega t}\frac{f(t)}{t-x}dt,\qquad \omega>0,\qquad x\geq 0,$$ where the bar indicates the Cauchy principal value and $f$ is a real-valued function with analytic continuation in the first quadrant, except possibly a branch point of algebraic type at the origin. When $x=0$, the integral is interpreted as a Hadamard finite-part integral, provided it is divergent. Asymptotic expansions in inverse powers of $\omega$ are derived for each fixed $x\geq 0$, which clarify the large $\omega$ behavior of this transform. We then present efficient and affordable approaches for numerical evaluation of such oscillatory transforms. Depending on the position of $x$, we classify our discussion into three regimes, namely, $x=\mathcal{O}(1)$ or $x\gg1$, $0<x\ll 1$ and $x=0$. Numerical experiments show that the convergence of the proposed methods greatly improve when the frequency $\omega$ increases. Some extensions to oscillatory Hilbert transforms with Bessel oscillators are briefly discussed as well.Comment: 32 pages, 6 figures, 4 table

Orthogonal polynomial expansions for the matrix exponential

AbstractMany different algorithms have been suggested for computing the matrix exponential. In this paper, we put forward the idea of expanding in either Chebyshev, Legendre or Laguerre orthogonal polynomials. In order for these expansions to converge quickly, we cluster the eigenvalues into diagonal blocks and accelerate using shifting and scaling

Spectral methods for partial differential equations

Origins of spectral methods, especially their relation to the Method of Weighted Residuals, are surveyed. Basic Fourier, Chebyshev, and Legendre spectral concepts are reviewed, and demonstrated through application to simple model problems. Both collocation and tau methods are considered. These techniques are then applied to a number of difficult, nonlinear problems of hyperbolic, parabolic, elliptic, and mixed type. Fluid dynamical applications are emphasized

Chebyshev Series Expansion of Inverse Polynomials

An inverse polynomial has a Chebyshev series expansion 1/\sum(j=0..k)b_j*T_j(x)=\sum'(n=0..oo) a_n*T_n(x) if the polynomial has no roots in [-1,1]. If the inverse polynomial is decomposed into partial fractions, the a_n are linear combinations of simple functions of the polynomial roots. If the first k of the coefficients a_n are known, the others become linear combinations of these with expansion coefficients derived recursively from the b_j's. On a closely related theme, finding a polynomial with minimum relative error towards a given f(x) is approximately equivalent to finding the b_j in f(x)/sum_(j=0..k)b_j*T_j(x)=1+sum_(n=k+1..oo) a_n*T_n(x), and may be handled with a Newton method providing the Chebyshev expansion of f(x) is known.Comment: LaTeX2e, 24 pages, 1 PostScript figure. More references. Corrected typos in (1.1), (3.4), (4.2), (A.5), (E.8) and (E.11

Basic Methods for Computing Special Functions

This paper gives an overview of methods for the numerical evaluation of special functions, that is, the functions that arise in many problems from mathematical physics, engineering, probability theory, and other applied sciences. We consider in detail a selection of basic methods which are frequently used in the numerical evaluation of special functions: converging and asymptotic series, including Chebyshev expansions, linear recurrence relations, and numerical quadrature. Several other methods are available and some of these will be discussed in less detail. We give examples of recent software for special functions where these methods are used. We mention a list of new publications on computational aspects of special functions available on our website

Hybrid Radial hp/Angular Galerkin Methods in Linearised Rotating Magnetohydrodynamics in Spheres and Related Geometries.

This thesis investigates hp-methods, involving refinement of mesh (h) and increase of polynomial order (p) used in several applications in linearised rotating magnetohydrodynamics in spherical geometries. The two hp-methods, one based on the Galerkin method, the other on a Chebyshev collocation method, are applied to eigenproblems with boundary layers for solving for viscous inertial modes and the onset of thermal convection at low Ekman number in a strong magnetic field, Elsasser number = 1. Results in one-dimension indicate the Chebyshev method is more efficient for forced problems, whereas the Galerkin method is more efficient for eigenproblems. Both hp-methods on a suitably defined radial mesh outperform the p-version for low Ekman number (10-6), and indicate increased efficiency as the Ekman number decreases. The Chebyshev method may find its use in nonlinear timestepping problems, which is beyond the scope of this thesis. A full sphere element is considered where it is important to build the analytic behaviour of the solutions at the origin into the approximation method. A method based on the weak Poincaré equation in a rotating tilted triaxial ellipsoid is presented. Extensions to ellipsoidal geometries are indicated