BZPEERPapers, people, evidence and research paths
Menu

Discover

SearchFind papers, people, topics, institutions and opportunities from one place.GraphSee how papers, authors, topics and evidence connect.

Research tools

HomeContinue from the papers, people and topics that matter now.CollectionsTurn saved papers into reading lists, evidence maps and shareable context.ResearchReconstruct the path that led you from search to paper to question.CollaborationMove from a paper, topic or expert request into a focused thread.PublishCreate full-text research articles with metadata, formulas and versions.

Network

ExpertsFind specialists whose work and reviews explain why they can help.CommunitiesJoin research groups organized around papers, topics and shared evidence.

Opportunities

ListingsBrowse roles, calls and collaborations connected to real research context.

Account

Log inReturn to your saved papers, graph views and active threads.Create accountStart saving papers, building a research trail and joining teams.
Log inCreate account

Source author record

Daan Huybrechs

Daan Huybrechs appears in the imported research catalog. Authorship, coauthor and topic links are available while profile ownership is still unclaimed.

ResearcherUnclaimed source record
math.NANumerical Analysismath.CAMathematical Software

Catalog footprint

What is connected

20works
4topics
4close collaborators

Actions

Connect this record

Log in to claim
Open graphBrowse works

Research graph

See the researcher in context

Open full explorer

Inspect adjacent papers, topics, institutions and collaborators without losing the researcher page.

Building this map preview

BZPEER is loading the nearby papers, people, topics and institutions for this page.

Published work

20 published item(s)

preprint2022arXiv

An analysis of least-squares oversampled collocation methods for compactly perturbed boundary integral equations in two dimensions

In recent work (Maierhofer & Huybrechs, 2022, Adv. Comput. Math.), the authors showed that least-squares oversampling can improve the convergence properties of collocation methods for boundary integral equations involving operators of certain pseudo-differential form. The underlying principle is that the discrete method approximates a Bubnov$-$Galerkin method in a suitable sense. In the present work, we extend this analysis to the case when the integral operator is perturbed by a compact operator $\mathcal{K}$ which is continuous as a map on Sobolev spaces on the boundary, $\mathcal{K}:H^{p}\rightarrow H^{q}$ for all $p,q\in\mathbb{R}$. This study is complicated by the fact that both the test and trial functions in the discrete Bubnov-Galerkin orthogonality conditions are modified over the unperturbed setting. Our analysis guarantees that previous results concerning optimal convergence rates and sufficient rates of oversampling are preserved in the more general case. Indeed, for the first time, this analysis provides a complete explanation of the advantages of least-squares oversampled collocation for boundary integral formulations of the Laplace equation on arbitrary smooth Jordan curves in 2D. Our theoretical results are shown to be in very good agreement with numerical experiments.

0 citations0 reviewsNo code linkedOpen access
preprint2022arXiv

Arbitrary-order asymptotic expansions of Gaussian quadrature rules with classical and generalised weight functions

Gaussian quadrature rules are a classical tool for the numerical approximation of integrals with smooth integrands and positive weight functions. We derive and expicitly list asymptotic expressions for the points and weights of Gaussian quadrature rules for three general classes of positive weight functions: analytic functions on a bounded interval with algebraic singularities at the endpoints, analytic weight functions on the halfline with exponential decay at infinity and an algebraic singularity at the finite endpoint, and analytic functions on the real line with exponential decay in both directions at infinity. The results include the Gaussian rules of classical orthogonal polynomials (Legendre, Jacobi, Laguerre and Hermite) as special cases. We present experiments indicating the range of the number of points at which these expressions achieve high precision. We provide an algorithm that can compute arbitrarily many terms in these expansions for the classical cases, and many though not all terms for the generalized cases.

0 citations0 reviewsNo code linkedOpen access
preprint2022arXiv

Frequency extraction for BEM-matrices arising from the 3D scalar Helmholtz equation

The discretisation of boundary integral equations for the scalar Helmholtz equation leads to large dense linear systems. Efficient boundary element methods (BEM), such as the fast multipole method (FMM) and $\Hmat$ based methods, focus on structured low-rank approximations of subblocks in these systems. It is known that the ranks of these subblocks increase linearly with the wavenumber. We explore a data-sparse representation of BEM-matrices valid for a range of frequencies, based on extracting the known phase of the Green's function. Algebraically, this leads to a Hadamard product of a frequency matrix with an $\Hmat$. We show that the frequency dependency of this $\Hmat$ can be determined using a small number of frequency samples, even for geometrically complex three-dimensional scattering obstacles. We describe an efficient construction of the representation by combining adaptive cross approximation with adaptive rational approximation in the continuous frequency dimension. We show that our data-sparse representation allows to efficiently sample the full BEM-matrix at any given frequency, and as such it may be useful as part of an efficient sweeping routine.

0 citations0 reviewsNo code linkedOpen access
preprint2022arXiv

On the computation of the SVD of Fourier submatrices

Contiguous submatrices of the Fourier matrix are known to be ill-conditioned. In a recent paper in SIAM Review A. Barnett has provided new bounds on the rate of ill-conditioning of the discrete Fourier submatrices. In this paper we focus on the corresponding singular value decomposition. The singular vectors go by the name of periodic discrete prolate spheroidal sequences (P-DPSS). The singular values exhibit an initial plateau, which depends on the dimensions of the submatrix, after which they decay rapidly. The latter regime is known as the plunge region and it is compatible with the submatrices being ill-conditioned. The discrete prolate sequences have received much less study than their continuous counterparts, prolate spheroidal wave functions, associated with continuous Fourier transforms and widely studied following the work of Slepian in the 1970's. In this paper we collect and expand known results on the stable numerical computation of the singular values and vectors of Fourier submatrices. We illustrate the computations and point out a few applications in which Fourier submatrices arise.

0 citations0 reviewsNo code linkedOpen access
preprint2020arXiv

Approximating smooth, multivariate functions on irregular domains

In this paper, we introduce a method known as polynomial frame approximation for approximating smooth, multivariate functions defined on irregular domains in $d$ dimensions, where $d$ can be arbitrary. This method is simple, and relies only on orthogonal polynomials on a bounding tensor-product domain. In particular, the domain of the function need not be known in advance. When restricted to a subdomain, an orthonormal basis is no longer a basis, but a frame. Numerical computations with frames present potential difficulties, due to the near-linear dependence of the truncated approximation system. Nevertheless, well-conditioned approximations can be obtained via regularization, for instance, truncated singular value decompositions. We comprehensively analyze such approximations in this paper, providing error estimates for functions with both classical and mixed Sobolev regularity, with the latter being particularly suitable for higher-dimensional problems. We also analyze the sample complexity of the approximation for sample points chosen randomly according to a probability measure, providing estimates in terms of the corresponding \textit{Nikolskii inequality} for the domain. In particular, we show that the sample complexity for points drawn from the uniform measure is quadratic (up to a log factor) in the dimension of the polynomial space, independently of $d$, for a large class of nontrivial domains. This extends a well-known result for polynomial approximation in hypercubes.

0 citations0 reviewsNo code linkedOpen access
preprint2020arXiv

Efficient function approximation on general bounded domains using splines on a Cartesian grid

Functions on a bounded domain in scientific computing are often approximated using piecewise polynomial approximations on meshes that adapt to the shape of the geometry. We study the problem of function approximation using splines on a regular but oversampled grid that is defined on a bounding box. This approach allows the use of high order and highly structured splines as a basis for piecewise polynomials. The methodology is analogous to that of Fourier extensions, using Fourier series on a bounding box, which leads to spectral accuracy for smooth functions. However, Fourier extension approximations involve solving a highly ill-conditioned linear system, and this is an expensive step. The computational complexity of recent algorithms is $\mathcal O\left(N\log^2\left(N\right)\right)$ in 1-D and $\mathcal O\left(N^2\log^2\left(N\right)\right)$ in 2-D. We show that, compared to Fourier extension, the compact support of B-splines enables improved complexity for multivariate approximations, namely $\mathcal O(N)$ in 1-D, $\mathcal O\left(N^{3/2}\right)$ in 2-D and more generally $\mathcal O\left(N^{3(d-1)/d}\right)$ in $d$-D with $d>1$. By using a direct sparse QR solver for a related linear system, we also observe that the computational complexity can be nearly linear in practice. This comes at the cost of achieving only algebraic rates of convergence. Our statements are corroborated with numerical experiments and Julia code is available.

0 citations0 reviewsNo code linkedOpen access
preprint2020arXiv

Efficient function approximation on general bounded domains using wavelets on a cartesian grid

Fourier extension is an approximation method that alleviates the periodicity requirements of Fourier series and avoids the Gibbs phenomenon when approximating functions. We describe a similar extension approach using regular wavelet bases on a hypercube to approximate functions on subsets of that cube. These subsets may have a general shape. This construction is inherently associated with redundancy which leads to severe ill-conditioning, but recent theory shows that nevertheless high accuracy and numerical stability can be achieved using regularization and oversampling. Regularized least squares solvers, such as the truncated singular value decomposition, that are suited to solve the resulting ill-conditioned and skinny linear system generally have cubic computational cost. We compare several algorithms that improve on this complexity. The improvements benefit from the sparsity in and the structure of the discrete wavelet transform. We present a method that requires $\mathcal O(N)$ operations in 1-D and $\mathcal O(N^{3(d-1)/d})$ in $d$-D, $d>1$. We experimentally show that direct sparse QR solvers appear to be more time-efficient, but yield larger expansion coefficients.

0 citations0 reviewsNo code linkedOpen access
preprint2020arXiv

Fast hybrid numerical-asymptotic boundary element methods for high frequency screen and aperture problems based on least-squares collocation

We present a hybrid numerical-asymptotic (HNA) boundary element method (BEM) for high frequency scattering by two-dimensional screens and apertures, whose computational cost to achieve any prescribed accuracy remains bounded with increasing frequency. Our method is a collocation implementation of the high order hp HNA approximation space of Hewett et al. IMA J. Numer. Anal. 35 (2015), pp.1698- 1728, where a Galerkin implementation was studied. An advantage of the current collocation scheme is that the one-dimensional highly oscillatory singular integrals appearing in the BEM matrix entries are significantly easier to evaluate than the two-dimensional integrals appearing in the Galerkin case, which leads to much faster computation times. Here we compute the required integrals at frequency-independent cost using the numerical method of steepest descent, which involves complex contour deformation. The change from Galerkin to collocation is nontrivial because naive collocation implementations based on square linear systems suffer from severe numerical instabilities associated with the numerical redundancy of the HNA basis, which produces highly ill-conditioned BEM matrices. In this paper we show how these instabilities can be removed by oversampling, and solving the resulting overdetermined collocation system in a weighted least-squares sense using a truncated singular value decomposition. On the basis of our numerical experiments, the amount of oversampling required to stabilise the method is modest (around 25% typically suffices) and independent of frequency. As an application of our method we present numerical results for high frequency scattering by prefractal approximations to the middle-third Cantor set.

0 citations0 reviewsNo code linkedOpen access
preprint2020arXiv

Frames and numerical approximation II: generalized sampling

In a previous paper [Adcock & Huybrechs, 2019] we described the numerical approximation of functions using redundant sets and frames. Redundancy in the function representation offers enormous flexibility compared to using a basis, but ill-conditioning often prevents the numerical computation of best approximations. We showed that, in spite of said ill-conditioning, approximations with regularization may still provide accuracy up to order $\sqrtε$, where $ε$ is a small truncation threshold. When using frames, i.e. complete systems that are generally redundant but which provide infinite representations with coefficients of bounded norm, this accuracy can actually be achieved for all functions in a space. Here, we generalize that setting in two ways. We assume information or samples from $f$ from a wide class of linear operators acting on $f$, rather than inner products associated with the best approximation projection. This enables the analysis of fully discrete approximations based, for instance, on function values only. Next, we allow oversampling, leading to least-squares approximations. We show that this leads to much improved accuracy on the order of $ε$ rather than $\sqrtε$. Overall, we demonstrate that numerical function approximation using redundant representations may lead to highly accurate approximations in spite of having to solve ill-conditioned systems of equations.

0 citations0 reviewsNo code linkedOpen access
preprint2020arXiv

On the adaptive spectral approximation of functions using redundant sets and frames

The approximation of smooth functions with a spectral basis typically leads to rapidly decaying coefficients where the rate of decay depends on the smoothness of the function and vice-versa. The optimal number of degrees of freedom in the approximation can be determined with relative ease by truncating the coefficients once a threshold is reached. Recent approximation schemes based on redundant sets and frames extend the applicability of spectral approximations to functions defined on irregular geometries and to certain non-smooth functions. However, due to their inherent redundancy, the expansion coefficients in frame approximations do not necessarily decay even for very smooth functions. In this paper, we highlight this lack of equivalence between smoothness and coefficient decay and we explore approaches to determine an optimal number of degrees of freedom for such redundant approximations.

0 citations0 reviewsNo code linkedOpen access
preprint2019arXiv

A numerical method for oscillatory integrals with coalescing saddle points

The value of a highly oscillatory integral is typically determined asymptotically by the behaviour of the integrand near a small number of critical points. These include the endpoints of the integration domain and the so-called stationary points or saddle points -- roots of the derivative of the phase of the integrand -- where the integrand is locally non-oscillatory. Modern methods for highly oscillatory quadrature exhibit numerical issues when two such saddle points coalesce. On the other hand, integrals with coalescing saddle points are a classical topic in asymptotic analysis, where they give rise to uniform asymptotic expansions in terms of the Airy function. In this paper we construct Gaussian quadrature rules that remain uniformly accurate when two saddle points coalesce. These rules are based on orthogonal polynomials in the complex plane. We analyze these polynomials, prove their existence for even degrees, and describe an accurate and efficient numerical scheme for the evaluation of oscillatory integrals with coalescing saddle points.

0 citations0 reviewsNo code linkedOpen access
preprint2015arXiv

Construction and implementation of asymptotic expansions for Jacobi--type orthogonal polynomials

We are interested in the asymptotic behavior of orthogonal polynomials of the generalized Jacobi type as their degree $n$ goes to $\infty$. These are defined on the interval $[-1,1]$ with weight function $w(x)=(1-x)^α(1+x)^βh(x)$, $α,β>-1$ and $h(x)$ a real, analytic and strictly positive function on $[-1,1]$. This information is available in the work of Kuijlaars, McLaughlin, Van Assche and Vanlessen, where the authors use the Riemann--Hilbert formulation and the Deift--Zhou non-linear steepest descent method. We show that computing higher-order terms can be simplified, leading to their efficient construction. The resulting asymptotic expansions in every region of the complex plane are implemented both symbolically and numerically, and the code is made publicly available. The main advantage of these expansions is that they lead to increasing accuracy for increasing degree of the polynomials, at a computational cost that is actually independent of the degree. In contrast, the typical use of the recurrence relation for orthogonal polynomials in computations leads to a cost that is at least linear in the degree. Furthermore, the expansions may be used to compute Gaussian quadrature rules in $\mathcal{O}(n)$ operations, rather than $\mathcal{O}(n^2)$ based on the recurrence relation.

0 citations0 reviewsNo code linkedOpen access
preprint2015arXiv

Fast Algorithms for the computation of Fourier Extensions of arbitrary length

Fourier series of smooth, non-periodic functions on $[-1,1]$ are known to exhibit the Gibbs phenomenon, and exhibit overall slow convergence. One way of overcoming these problems is by using a Fourier series on a larger domain, say $[-T,T]$ with $T>1$, a technique called Fourier extension or Fourier continuation. When constructed as the discrete least squares minimizer in equidistant points, the Fourier extension has been shown shown to converge geometrically in the truncation parameter $N$. A fast ${\mathcal O}(N \log^2 N)$ algorithm has been described to compute Fourier extensions for the case where $T=2$, compared to ${\mathcal O}(N^3)$ for solving the dense discrete least squares problem. We present two ${\mathcal O}(N\log^2 N )$ algorithms for the computation of these approximations for the case of general $T$, made possible by exploiting the connection between Fourier extensions and Prolate Spheroidal Wave theory. The first algorithm is based on the explicit computation of so-called periodic discrete prolate spheroidal sequences, while the second algorithm is purely algebraic and only implicitly based on the theory.

0 citations0 reviewsNo code linkedOpen access
preprint2015arXiv

Fast and highly accurate computation of Chebyshev expansion coefficients of analytic functions

Chebyshev expansion coefficients can be computed efficiently by using the FFT, and for smooth functions the resulting approximation is close to optimal, with computations that are numerically stable. Given sufficiently accurate function samples, the Chebyshev expansion coefficients can be computed to machine precision accuracy. However, the accuracy is only with respect to absolute error, and this implies that very small expansion coefficients typically have very large relative error. Upon differentiating a Chebyshev expansion, this relative error in the small coefficients is magnified and accuracy may be lost, especially after repeated differentiation. At first sight, this seems unavoidable. Yet, in this paper, we focus on an alternative computation of Chebyshev expansion coefficients using contour integrals in the complex plane. The main result is that the coefficients can be computed with machine precision relative error, rather than absolute error. This implies that even very small coefficients can be computed with full floating point accuracy, even when they are themselves much smaller than machine precision. As a result, no accuracy is lost after differentiating the expansion, and even the 100th derivative of an analytic function can be computed with near machine precision accuracy using standard floating point arithmetic. In some cases, the contour integrals can be evaluated using the FFT, making the approach both highly accurate and fast.

0 citations0 reviewsNo code linkedOpen access
preprint2013arXiv

On the numerical stability of Fourier extensions

An effective means to approximate an analytic, nonperiodic function on a bounded interval is by using a Fourier series on a larger domain. When constructed appropriately, this so-called Fourier extension is known to converge geometrically fast in the truncation parameter. Unfortunately, computing a Fourier extension requires solving an ill-conditioned linear system, and hence one might expect such rapid convergence to be destroyed when carrying out computations in finite precision. The purpose of this paper is to show that this is not the case. Specifically, we show that Fourier extensions are actually numerically stable when implemented in finite arithmetic, and achieve a convergence rate that is at least superalgebraic. Thus, in this instance, ill-conditioning of the linear system does not prohibit a good approximation. In the second part of this paper we consider the issue of computing Fourier extensions from equispaced data. A result of Platte, Trefethen & Kuijlaars states that no method for this problem can be both numerically stable and exponentially convergent. We explain how Fourier extensions relate to this theoretical barrier, and demonstrate that they are particularly well suited for this problem: namely, they obtain at least superalgebraic convergence in a numerically stable manner.

0 citations0 reviewsNo code linkedOpen access
preprint2013arXiv

Zero distribution of complex orthogonal polynomials with respect to exponential weights

We study the limiting zero distribution of orthogonal polynomials with respect to some particular exponential weights exp(-nV(z)) along contours in the complex plane. We are especially interested in the question under which circumstances the zeros of the orthogonal polynomials accumulate on a single analytic arc (one cut case), and in which cases they do not. In a family of cubic polynomial potentials V(z) = - iz^3/3 + iKz, we determine the precise values of K for which we have the one cut case. We also prove the one cut case for a monomial quintic V(z) = - iz^5/5 on a contour that is symmetric in the imaginary axis.

0 citations0 reviewsNo code linkedOpen access
preprint2012arXiv

A Gaussian quadrature rule for oscillatory integrals on a bounded interval

We investigate a Gaussian quadrature rule and the corresponding orthogonal polynomials for the oscillatory weight function $e^{iωx}$ on the interval $[-1,1]$. We show that such a rule attains high asymptotic order, in the sense that the quadrature error quickly decreases as a function of the frequency $ω$. However, accuracy is maintained for all values of $ω$ and in particular the rule elegantly reduces to the classical Gauss-Legendre rule as $ω\to 0$. The construction of such rules is briefly discussed, and though not all orthogonal polynomials exist, it is demonstrated numerically that rules with an even number of points are always well defined. We show that these rules are optimal both in terms of asymptotic order as well as in terms of polynomial order.

0 citations0 reviewsNo code linkedOpen access
preprint2012arXiv

Explicit barycentric weights for polynomial interpolation in the roots or extrema of classical orthogonal polynomials

Barycentric interpolation is arguably the method of choice for numerical polynomial interpolation. The polynomial interpolant is expressed in terms of function values using the so-called barycentric weights, which depend on the interpolation points. Few explicit formulae for these barycentric weights are known. In [H. Wang and S. Xiang, Math. Comp., 81 (2012), 861--877], the authors have shown that the barycentric weights of the roots of Legendre polynomials can be expressed explicitly in terms of the weights of the corresponding Gaussian quadrature rule. This idea was subsequently implemented in the Chebfun package [L. N. Trefethen and others, The Chebfun Development Team, 2011] and in the process generalized by the Chebfun authors to the roots of Jacobi, Laguerre and Hermite polynomials. In this paper, we explore the generality of the link between barycentric weights and Gaussian quadrature and show that such relationships are related to the existence of lowering operators for orthogonal polynomials. We supply an exhaustive list of cases, in which all known formulae are recovered and also some new formulae are derived, including the barycentric weights for Gauss-Radau and Gauss-Lobatto points. Based on a fast ${\mathcal O}(n)$ algorithm for the computation of Gaussian quadrature, due to Hale and Townsend, this leads to an ${\mathcal O}(n)$ computational scheme for barycentric weights.

0 citations0 reviewsNo code linkedOpen access
preprint2012arXiv

On the resolution power of Fourier extensions for oscillatory functions

Functions that are smooth but non-periodic on a certain interval possess Fourier series that lack uniform convergence and suffer from the Gibbs phenomenon. However, they can be represented accurately by a Fourier series that is periodic on a larger interval. This is commonly called a Fourier extension. When constructed in a particular manner, Fourier extensions share many of the same features of a standard Fourier series. In particular, one can compute Fourier extensions which converge spectrally fast whenever the function is smooth, and exponentially fast if the function is analytic, much the same as the Fourier series of a smooth/analytic and periodic function. With this in mind, the purpose of this paper is to describe, analyze and explain the observation that Fourier extensions, much like classical Fourier series, also have excellent resolution properties for representing oscillatory functions. The resolution power, or required number of degrees of freedom per wavelength, depends on a user-controlled parameter and, as we show, it varies between 2 and π. The former value is optimal and is achieved by classical Fourier series for periodic functions, for example. The latter value is the resolution power of algebraic polynomial approximations. Thus, Fourier extensions with an appropriate choice of parameter are eminently suitable for problems with moderate to high degrees of oscillation.

0 citations0 reviewsNo code linkedOpen access
preprint2011arXiv

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ωt})(x)=-int_{0}^{\infty}e^{iωt}\frac{f(t)}{t-x}dt,\qquad ω>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 $ω$ are derived for each fixed $x\geq 0$, which clarify the large $ω$ 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 $ω$ increases. Some extensions to oscillatory Hilbert transforms with Bessel oscillators are briefly discussed as well.

0 citations0 reviewsNo code linkedOpen access