A Riemann--Hilbert approach to computing the inverse spectral map for measures supported on disjoint intervals

TL;DR

This paper introduces a fast numerical method based on Riemann--Hilbert problems for computing orthogonal polynomials on disjoint intervals, enabling efficient calculation of recurrence coefficients and polynomial evaluations.

Contribution

It develops an $ ext{O}(N)$ Riemann--Hilbert based approach for orthogonal polynomials on multiple intervals, addressing a gap where analytical formulas are unknown.

Findings

Method outperforms existing techniques in efficiency.

Demonstrates fast convergence in applications.

Enables pointwise evaluation of polynomials and Cauchy transforms.

Abstract

We develop a numerical method for computing with orthogonal polynomials that are orthogonal on multiple, disjoint intervals for which analytical formulae are currently unknown. Our approach exploits the Fokas--Its--Kitaev Riemann--Hilbert representation of the orthogonal polynomials to produce an $\mathrm{O}(N)$ method to compute the first $N$ recurrence coefficients. The method can also be used for pointwise evaluation of the polynomials and their Cauchy transforms throughout the complex plane. The method encodes the singularity behavior of weight functions using weighted Cauchy integrals of Chebyshev polynomials. This greatly improves the efficiency of the method, outperforming other available techniques. We demonstrate the fast convergence of our method and present applications to integrable systems and approximation theory.

Equations

$$w(x)=\sum _{j=1}^{g+1}\text{mathbbm}{1}_{[a_j,b_j]}(x)h_j(x){\left(\sqrt{x-a_j}\right)}^{{\alpha }_j}{\left(\sqrt{b_j-x}\right)}^{{\beta }_j},{\alpha }_j,{\beta }_j\in \{-1,1\},$$

$${\int }_{\Sigma }{\pi }_j(x){\pi }_k(x)w(x)\mathrm{d}x=h_j{\delta }_{jk},$$

$$p_j(x)={\gamma }_j{\pi }_j(x),$$

$${\int }_{\Sigma }{p}_j(x)p_k(x)w(x)\mathrm{d}x={\delta }_{jk}.$$

$$x{p}_0(x)=a_0{p}_0(x)+b_0{p}_1(x),$$

$$x{p}_k(x)=b_{k-1}{p}_{k-1}(x)+a_k{p}_k(x)+b_k{p}_{k+1}(x),k\ge 1,$$

$$\mathbf{J}(w)=\left(\begin{array}{ccccc}{a}_0& b_0& & & \\ b_0& a_1& b_1& & \\ & b_1& a_2& b_2& \\ & & \ddots & \ddots & \ddots \end{array}\right).$$

$$\mathbf{p}(x)=\left(\begin{array}{cccc}{p}_0(x)& p_1(x)& p_2(x)& \cdots \end{array}\right),$$

$$\mathbf{p}(x{)}^{T}x=\mathbf{J}(w)\mathbf{p}(x{)}^T.$$

$$w_{\mathcal{T}}(x)=\frac{1}{\pi \sqrt{1-x}\sqrt{x+1}},$$

$$\mathbf{J}(w_{\mathcal{T}})=\left(\begin{array}{ccccc}0& 1/\sqrt{2}& & & \\ 1/\sqrt{2}& 0& 1/2& & \\ & 1/2& 0& 1/2& \\ & & \ddots & \ddots & \ddots \end{array}\right).$$

$$w_U(x)=\frac{2}{\pi }\sqrt{1-x}\sqrt{x+1},$$

$$\mathbf{J}(w_U)=\left(\begin{array}{ccccc}0& 1/2& & & \\ 1/2& 0& 1/2& & \\ & 1/2& 0& 1/2& \\ & & \ddots & \ddots & \ddots \end{array}\right).$$

$$w_V(x)=\frac{1}{\pi }\frac{\sqrt{x+1}}{\sqrt{1-x}},$$

$$\mathbf{J}(w_V)=\left(\begin{array}{ccccc}1/2& 1/2& & & \\ 1/2& 0& 1/2& & \\ & 1/2& 0& 1/2& \\ & & \ddots & \ddots & \ddots \end{array}\right).$$

$$w_W(x)=\frac{1}{\pi }\frac{\sqrt{1-x}}{\sqrt{x+1}},$$

$$\mathbf{J}(w_W)=\left(\begin{array}{ccccc}-1/2& 1/2& & & \\ 1/2& 0& 1/2& & \\ & 1/2& 0& 1/2& \\ & & \ddots & \ddots & \ddots \end{array}\right).$$

$${\mathcal{C}}_{\Gamma }f(z)=\frac{1}{2\pi \mathrm{i}}{\int }_{\Gamma }\frac{f(s)}{s-z}\mathrm{d}s.$$

$$F^{\pm }(z)=\underset{\begin{array}{c}{z}^{\prime }\to z,\\ z\in B_{\pm }\end{array}}{\lim }F(z^{\prime }),$$

$${\mathcal{C}}_{\Gamma }^{\pm }f(z)=\underset{\begin{array}{c}{z}^{\prime }\to z,\\ z\in B_{\pm }\end{array}}{\lim }{\mathcal{C}}_{\Gamma }f(z^{\prime }).$$

$$C_k(z)={\mathcal{C}}_{\Sigma }[p_{k}w](z)=\frac{1}{2\pi \mathrm{i}}{\int }_{\Sigma }\frac{p_k(s)w(s)}{s-z}\mathrm{d}s,$$

$$z{C}_0(z)=a_0{C}_0(z)+b_0{C}_1(z)-\frac{1}{2\pi \mathrm{i}}{\int }_{\Sigma }w(x)\mathrm{d}x,$$

$$z{C}_k(z)=b_{k-1}{C}_{k-1}(z)+a_k{C}_k(z)+b_k{C}_{k+1}(z),k\ge 1.$$

$${\mathbf{Y}}_n\text{leavevmode}\text{is analytic in}\text{leavevmode}\mathbb{C}\setminus ,$$

$${\mathbf{Y}}_n^{+}(z)={\mathbf{Y}}_n^{-}(z)\left(\begin{array}{cc}1& w(z)\\ 0& 1\end{array}\right),z\in ,$$

$${\mathbf{Y}}_n(z)=\left(\mathbf{I}+O(z^{-1})\right)\left(\begin{array}{cc}{z}^n& 0\\ 0& z^{-n}\end{array}\right),z\to \infty .$$

$${\mathbf{Y}}_n(z)=\left(\begin{array}{cc}{\pi }_n(z)& {\mathcal{C}}_{\Sigma }[{\pi }_{n}w](z)\\ -2\pi \mathrm{i}{\gamma }_{n-1}^2{\pi }_{n-1}(z)& -2\pi \mathrm{i}{\gamma }_{n-1}^2{\mathcal{C}}_{\Sigma }[{\pi }_{n-1}w](z)\end{array}\right).$$

$${\mathbf{Y}}_n\text{leavevmode}\text{is analytic in}\text{leavevmode}\mathbb{C}\setminus \Sigma ,$$

$${\mathbf{Y}}_n^{+}(z)={\mathbf{Y}}_n^{-}(z)\left(\begin{array}{cc}1& w(z)\\ 0& 1\end{array}\right),z\in \Sigma ,$$

$${\mathbf{Y}}_n(z)=\left(\mathbf{I}+O(z^{-1})\right)\left(\begin{array}{cc}{z}^n& 0\\ 0& z^{-n}\end{array}\right),z\to \infty ,$$

$${\mathbf{Y}}_n(z)=O\left(\begin{array}{cc}1& 1+|z-a_j{|}^{{\alpha }_j/2}\\ 1& 1+|z-a_j{|}^{{\alpha }_j/2}\end{array}\right),z\to a_j,$$

$${\mathbf{Y}}_n(z)=O\left(\begin{array}{cc}1& 1+|z-b_j{|}^{{\beta }_j/2}\\ 1& 1+|z-b_j{|}^{{\beta }_j/2}\end{array}\right),z\to b_j,$$

$$w_j(x)=h_j(x){\left(\sqrt{x-a_j}\right)}^{{\alpha }_j}{\left(\sqrt{b_j-x}\right)}^{{\beta }_j}.$$

$${\mathbf{Z}}_n(z)=\{\begin{array}{ll}{\mathbf{Y}}_n(z)\left(\begin{array}{cc}1& 0\\ -1/w_j(z)& 1\end{array}\right)& z\in D_j\cap {\mathbb{C}}^{+},\\ {\mathbf{Y}}_n(z)\left(\begin{array}{cc}1& 0\\ 1/w_j(z)& 1\end{array}\right)& z\in D_j\cap {\mathbb{C}}^{-},\\ {\mathbf{Y}}_n(z)& \text{otherwise},\end{array}$$

Full text

This material is based on work supported by the National Science Foundation under Grant No. DMS-1945652 (TT). Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.

1. Introduction

Orthogonal polynomials are fundamental objects of study in computational mathematics with numerous applications ranging from quadrature [18] and eigenvalue approximation [26] to differential equations [32] and random matrix theory [9], and beyond [10, 15, 31, 36]. The theory of orthogonal polynomials dates back to the work of Chebyshev, Jacobi, Lagrange, and others. Much of this theory is summarized in the book of Szegő [36]. The modern theory of orthogonal polynomials has been intertwined with the theory of Riemann–Hilbert problems since the discovery of the Fokas, Its, and Kitaev (FIK) Riemann–Hilbert problem [14]. This characterization of orthogonal polynomials in terms of a solution of a Riemann–Hilbert problem, combined with the Deift–Zhou method of nonlinear steepest descent [11], has significantly advanced the asymptotic knowledge of orthogonal polynomials. Much of this advancement is presented in the book of Deift [9] which is fundamental to our work (see also [24]).

Orthogonal polynomials are arguably best utilized via the coefficients that make up their three-term recurrence. These recurrence coefficients, in a sense, fully characterize their orthogonal polynomials, allowing one to compute their values, their zeroes, and their normalization coefficients [15]. Thus, orthogonal polynomials that are otherwise not characterized may be utilized numerically if their recurrence coefficients can be reliably computed. The Jacobi matrix of three-term recurrence coefficients has the orthogonality measure as its spectral measure [9], and so computing the three-term recurrence coefficients, given the measure, is often referred to as the computing the inverse spectral map.

Historically, given an orthogonality weight, the orthogonal polynomial recurrence coefficients have been computed via a discretization of the Gram–Schmidt procedure known as the Stieltjes procedure [15]. While reliable, this method, in a different setting, has been realized to require $\operatorname{O}(N^{3})$ arithmetic operations to compute the first $N$ recurrence coefficients [38], leaving significant room for speedup through a Riemann–Hilbert approach. An $\operatorname{O}(N^{2})$ algorithm based on [19] is presented in [15] but is not without its own limitations. Obtaining numerical solutions to the FIK Riemann–Hilbert problem is a relatively recent area of research, and such methods are primarily based on techniques developed by Olver [30]. These techniques were applied in [38] to compute recurrence coefficients for orthogonal polynomials with exponential weights. This work demonstrated that a Riemann–Hilbert approach requires only $\operatorname{O}(N)$ arithmetic operations to compute the first $N$ recurrence coefficients. Furthermore, this approach allows such coefficients to be computed individually in $\operatorname{O}(1)$ arithmetic operations so that methods based on this approach are optimal in a time complexity sense. It is also possible that an increasingly analytical approach could be employed [23] to accomplish our same goal, and we leave this for future work, but the complication here is that asymptotic formulae are given in terms of theta functions [12, 40] which are themselves non-trivial to compute [8].

In this paper, we develop and apply a numerical Riemann–Hilbert approach to the problem presented in [25] and further extended in [12] to compute the recurrence coefficients for Chebyshev-like polynomials with support on multiple, disjoint intervals on the real line. Such orthogonal polynomials have both analytical and computational applications [17, 2, 34, 7]. For example, such orthogonal polynomials produce a new class of iterative solvers [34], improve matrix spectrum approximation when bounds on the eigenvalues are known a priori [6], and describe the numerical evolution of a Toda lattice [10]. Analytical formulae are not known for most of the cases we consider. The method remains accurate as the degree of the polynomial under consideration grows, and it bypasses the need to compute theta functions. More specifically, we compute orthogonal polynomial recurrence coefficients for weight functions of the form

[TABLE]

for a set of disjoint intervals, $\bigcup_{j=1}^{g+1}[a_{j},b_{j}]$ where each $h_{j}$ is positive on $[a_{j},b_{j}]$ and analytic in a neighborhood $\Omega_{j}$ of $[a_{j},b_{j}]$. We include a plot of a weight function of this form in Figure 1.

In most cases, it will suffice to consider $h_{j}(x)=1$ for all $j$, and recent work [20] has demonstrated that recurrence coefficients corresponding to a different choice of $h_{j}$ can be obtained from those corresponding to $h_{j}(x)=1$. In principle, our weight function (1) could be generalized to a Jacobi-like weight by allowing for any $\alpha_{j},\beta_{j}>-2$; however this would complicate our numerics significantly, as it would result in a Riemann–Hilbert problem defined on self-intersecting contours. Furthermore, the asymptotics are different if $\alpha_{j},\beta_{j}$ are not integers, with error terms behaving polynomially rather than exponentially [25]. Thus, while our method could, in principle, be extended to an $\operatorname{O}(N)$ method for Jacobi-like weights, this method would likely have a much larger leading constant than the method that we consider.

In Section 2, we introduce notation for orthogonal polynomials and Cauchy integrals. In Section 3, we apply a portion of the method of nonlinear steepest descent to the FIK Riemann–Hilbert problem with weight (1) and introduce an additional transformation so that a numerical solution can be computed. In Section 4, we develop numerical methods to solve the FIK Riemann–Hilbert problem and recover recurrence coefficients. These methods encode the singularity behavior of the weight (1), enabling much faster computations than the methods found in [38]. For example, our method computes the first 51 recurrence coefficients for a Chebyshev-$\mathcal{T}$-like weight ($\alpha_{j}=\beta_{j}=-1$) on the domain $[-1.8,-1]\cup[2,3]$ to machine precision in 0.47 seconds111All computations in this paper are performed on a Lenovo laptop running Ubuntu version 20.04 with 8 cores and 16 GB of RAM with an Intel® Core™ i7-11800H processor running at 2.30 GHz.. A further timing comparison with an optimized $\operatorname{O}(N^{2})$ method [19] can be found in Section 5.2. In Section 5, we demonstrate the numerical accuracy of our method and apply our method to evolving a Toda lattice and approximating $1/x$ on a disconnected domain. Code used to generate the plots in this paper can be found at [3].

Our method could also potentially be extended to compute the associated $N$-point Gauss quadrature nodes and weights for (1) in $\operatorname{O}(N)$ arithmetic operations using techniques similar to [21, 37, 5, 23]. Many extra complexities are introduced for (1), including the possibility of orthogonal polynomials having roots outside the support of the weight. For this approach to have any advantage over using the individual rules for each subinterval $[a_{j},b_{j}]$, it is likely that $g$ will need to be large.

2. Orthogonal polynomials and Cauchy integrals

2.1. Orthogonal polynomials

For our purposes, a weight function $w$ is a nonnegative function defined on $\Sigma\subset$ that is continuous and positive on the interior of $\Sigma$. Consider a sequence of univariate monic polynomials $\pi_{0}(x),\pi_{1}(x),\pi_{2}(x),\ldots$ such that $\pi_{j}$ has degree $j$ for all $j\in\mathbb{N}$. These polynomials are said to be orthogonal with respect to a weight function $w$ if

[TABLE]

where $h_{j}>0$ and $\delta_{jk}$ is the Kronecker delta.

The orthonormal polynomials $p_{0}(x),p_{1}(x),p_{2}(x),\ldots$ are defined such that

[TABLE]

for all $j\in\mathbb{N}$ where $\gamma_{j}=\frac{1}{\sqrt{h_{j}}}$. Clearly, the orthonormal polynomials satisfy

[TABLE]

Additionally, the orthonormal polynomials satisfy a symmetric three-term recurrence

[TABLE]

where $b_{k}>0$ for all $k$. The Jacobi operator associated with a sequence of orthonormal polynomials is defined as the infinite tridiagonal matrix

[TABLE]

If we let

[TABLE]

the three-term recurrence can be expressed as

[TABLE]

Chebyshev polynomials defined on $\mathbb{I}=[-1,1]$ will be particularly useful. A general reference is [27]. Normalized Chebyshev polynomials of the first kind are denoted by $\mathcal{T}_{k}(x)$ and have weight function

[TABLE]

and Jacobi operator

[TABLE]

We will always use the square root symbol to denote the principal branch defined on $\mathbb{C}\setminus(-\infty,0)$. Thus, $w_{\mathcal{T}}$ is analytic on $\mathbb{C}\setminus\left(\setminus(-1,1)\right).$ Note that our convention relates to the classical Chebyshev-$T$ polynomials by $\mathcal{T}_{0}(x)=T_{0}(x)$, and $\mathcal{T}_{j}(x)=\sqrt{2}T_{j}(x)$ for $j>0$.

Chebyshev polynomials of the second kind on $\mathbb{I}$ are denoted by $U_{k}(x)$ and have weight function

[TABLE]

and Jacobi operator

[TABLE]

Chebyshev polynomials of the third kind on $\mathbb{I}$ are denoted by $V_{k}(x)$ and have weight function

[TABLE]

and Jacobi operator

[TABLE]

Chebyshev polynomials of the fourth kind on $\mathbb{I}$ are denoted by $W_{k}(x)$ and have weight function

[TABLE]

and Jacobi operator

[TABLE]

2.2. Cauchy integrals

Given a curve $\Gamma\subset\mathbb{C}$ and a function $f:\Gamma\to\mathbb{C}$, the Cauchy transform $\mathcal{C}_{\Gamma}$ is as an operator that maps $f$ to its Cauchy integral, i.e.,

[TABLE]

If $\Gamma$ is a smooth, non-self-intersecting, oriented contour, then to each interior point $z\in\Gamma$, we associate a ball $B=B(z)$ such that $B\setminus\Gamma$ has two connected components $B_{\pm}$ where $B_{+}$ ($B_{-}$) is to the left (right) with respect to the orientation of $\Gamma$. For $F:\mathbb{C}\setminus\Gamma\to\mathbb{C}$, define the (non-tangential) boundary values by

[TABLE]

where this limit exists. We define the notation for the Cauchy operators

[TABLE]

Cauchy integrals of orthogonal polynomials will also be of use in this work. For orthonormal polynomials $p_{k}(x)$ associated with weight $w$ on $\Sigma\subset$, define

[TABLE]

for $z\in\mathbb{C}\setminus\Sigma$. It is well-known that if the three-term recurrence for $p_{k}$ is given by (2), these Cauchy integrals also satisfy the three-term recurrence [31]:

[TABLE]

3. The Riemann–Hilbert problem

3.1. Setting up the Riemann–Hilbert Problem

The original formulation of orthogonal polynomials as a matrix Riemann–Hilbert problem is due to Fokas, Its, and Kitaev and is for orthogonal polynomials with $\Sigma=$ [14]. Let $\pi_{0},\pi_{1},\ldots$ be the sequence of monic orthogonal polynomials for weight function $w$. Consider the Riemann–Hilbert problem for the unknown $\mathbf{Y}_{n}:\mathbb{C}\setminus\to\mathbb{C}^{2\times 2}$

[TABLE]

This problem has a unique solution that is given by [9]

[TABLE]

Then, (4) can be modified for orthogonal polynomials defined on a bounded set $\Sigma\subset$ by including asymptotic conditions at the boundary of $\Sigma$ [25]. We specifically consider orthogonal polynomials for a weight function (1) as discussed above. Recall that $\alpha_{j},\beta_{j}$ determine which Chebyshev weight function variant $w$ behaves like on $[a_{j},b_{j}]$.

The Riemann–Hilbert problem for this weight function is given by

[TABLE]

where the two asymptotic statements are taken element-wise. The unique solution is still given by (5). For computational ease, we do not require the weight function $w$ to be normalized.

3.2. Transforming the Riemann–Hilbert problem

We transform this problem into one that can be solved numerically through the following procedure:

• •

Step 1: Lens the Riemann–Hilbert problem, factoring and moving jump conditions to regions where we will induce exponential decay.

• •

Step 2: Determine and apply a differential to remove singularities at infinity.

• •

Step 3: Determine and apply an auxiliary function to remove the jump conditions on the gaps induced by the differential, and “solve” the resulting problem.

This procedure is performed in full for orthogonal polynomials with support on a single interval in [25, 24] and for multiple intervals in [12, 40]; the novelty of our approach lies in how Step 3 is applied and the fact that “solve” means constructing a reliable numerical approximation, not an asymptotic approximation.

3.2.1. Step 1: Lensing

Consider

[TABLE]

This function is analytic on $\Omega_{j}\setminus(\setminus(a_{j},b_{j}))$. Let $C_{j}$ be a counterclockwise-oriented curve lying in $\Omega_{j}$ that encircles, but does not intersect $[a_{j},b_{j}]$, and denote its interior by $D_{j}$. In the end, our numerical method will take $C_{j}$ to be a circle, imposing additional restrictions on $\Omega_{j}$. We include examples of these regions in Figure 2.

Define the transformation

[TABLE]

for $j=1,\ldots,g+1$. It follows that for $z\in C_{j}\cap\mathbb{C}^{+}$,

[TABLE]

and for $z\in C_{j}\cap\mathbb{C}^{-}$,

[TABLE]

For $z\in(a_{j},b_{j})$,

[TABLE]

Now, we compute the jump condition on $D_{j}\cap(\setminus(a_{j},b_{j}))$. For $x>b_{j}$ or $x<a_{j}$, $w_{j}^{+}(x)=-w_{j}^{-}(x)$ where now the $+$ denotes a limit from the upper half plane, and $-$ denotes a limit from the lower. Then,

[TABLE]

We find that there is no additional jump condition induced between $[a_{j},b_{j}]$ and $C_{j}$. This is a special situation in which no local parametrices are required to compute the asymptotics of the polynomials. The asymptotic behavior at the endpoints is multiplied by $\pm 1/w_{j}(z)$ in this transformation, and the lensed Riemann–Hilbert problem is given by

[TABLE]

3.2.2. Step 2: The differential $\mathfrak{g}^{\prime}$

To remove the growth/decay at infinity, we seek a function $\mathfrak{g}$ that satisfies the following:

• •

$\mathfrak{g}^{\prime}(z)=z^{-1}+\operatorname{O}(z^{-2})$ as $z\to\infty$.

• •

$\mathfrak{g}^{\prime+}(z),\mathfrak{g}^{\prime-}(z)$ are purely imaginary on $[a_{j},b_{j}]$ for $j=1,\ldots,g+1$.

• •

$\int_{b_{j}}^{a_{j+1}}\mathfrak{g}^{\prime}(z)\mathrm{d}z=0$ for $j=1,\ldots,g$.

This is referred to as the exterior Green’s function with pole at infinity. We achieve this by defining

[TABLE]

Here, $Q_{g}$ is a monic polynomial of degree $g$

[TABLE]

whose coefficients $h_{k}$ are determined by the linear system

[TABLE]

for $j=1,\ldots,g$. The unique solvability of this system follows from the fact that $\frac{z^{k}}{R(z)}$ for $k=0,1,\ldots,g-1$ form a basis for holomorphic differentials on the hyperelliptic Riemann surface defined by $w^{2}=R(z)^{2}$. Indeed, the coefficient matrix here is $-\frac{1}{2}\tilde{\mathbf{A}}$ where $\tilde{\mathbf{A}}$ is invertible as in Lemma 3.1 below. Then, we define $\mathfrak{g}$ by normalizing at the leftmost endpoint. If some function satisfies $\psi^{\prime}=\mathfrak{g}^{\prime}$, then we let

[TABLE]

We also define

[TABLE]

for $z\in(b_{j},a_{j+1})$, $j=1,\ldots,g$, where this value is independent of $z$. We have that

[TABLE]

for $z\in(a_{j},b_{j})$. Now, we consider the transformation

[TABLE]

where $\sigma_{3}$ is the third Pauli matrix, i.e., $\mathfrak{g}$ has asymptotic expansion

[TABLE]

as $z\to\infty$. This asymptotic form follows from the asymptotic expansion of $\mathfrak{g}^{\prime}$ at infinity. By construction, $\mathbf{S}_{n}(z)=\mathbf{I}+\operatorname{O}(z^{-1})$ as $z\to\infty$. Now, the jump condition for $z\in C_{j}\cap\mathbb{C}^{+}$ is given by

[TABLE]

and similarly, the jump condition for $z\in C_{j}\cap\mathbb{C}^{-}$ is given by

[TABLE]

For $z\in(a_{j},b_{j})$, the jump condition is unchanged:

[TABLE]

For $z\in(b_{j},a_{j+1})$, the jump condition is given by

[TABLE]

The asymptotic conditions at the endpoints remain unchanged since $|\mathfrak{g}|$ is bounded in the neighborhood of the contours. We have arrived at the normalized Riemann–Hilbert problem:

[TABLE]

Note that $\operatorname{Re}\mathfrak{g}(z)>0$ for $z\in C_{j}$ (see [12, Appendix B], for example) so the jump condition on $C_{j}$ tends exponentially to the identity as $n\to\infty$.

3.2.3. Step 3: Removing the jump conditions on the gaps

Since the normalized Riemann–Hilbert problem (11) contains jump conditions on the gaps between intervals, the jump matrix is discontinuous. In order to appropriately solve this problem numerically, we must further transform it. We do this by devising an auxiliary function to remove the jump conditions on the gaps between intervals.

Consider a function $\mathfrak{h}_{n}$ defined as

[TABLE]

where $A_{1}(n),\ldots,A_{g+1}(n)$ are undetermined constants. Recall that on the gaps $(b_{j},a_{j+1})$, $R$ is analytic. Since $R(z)\sim z^{g+1}$ as $z\to\infty$, we can use the asymptotics of Cauchy integrals to obtain conditions on $A_{1}(n),\ldots,A_{g+1}(n)$ such that $\mathfrak{h}_{n}(z)\sim z^{-1}$ as $z\to\infty$. We need to solve

[TABLE]

for $k=1,\ldots,g+1$ for constants $A_{1}(n),\ldots,A_{g+1}(n)$. This forms an $(g+1)\times(g+1)$ linear system. Note that $A_{1}(n),\ldots,A_{g+1}(n)\in$ and are bounded above and below due to the fact that $\Delta_{j}$ is purely imaginary, and the coefficient matrix does not depend on $n$. Furthermore, this coefficient matrix is always invertible.

Lemma 3.1**.**

The $(g+1)\times(g+1)$ matrix $\mathbf{H}$ defined by

[TABLE]

is invertible.

Proof.

As noted before, $\frac{z^{k-1}}{R(z)}$ for $k=1,\ldots,g$ form a basis for holomorphic differentials on the hyperelliptic Riemann surface defined by $w^{2}=R(z)^{2}$. Furthermore, for a differential $\omega$, if we define the cycle $\tilde{\mathfrak{b}}_{j}$ by

[TABLE]

the classical $b$-periods, denoted by $\mathfrak{b}_{j}$, $j=1,2,\ldots,g$, can be expressed via

[TABLE]

Define the $a$- and $b$-period matrices $\tilde{\mathbf{A}}$ and $\tilde{\mathbf{B}}$, respectively, by

[TABLE]

Standard theory implies that $\tilde{\mathbf{A}}$ is invertible. This implies that if $\mathbf{H}_{g}$ is the upper-left $g\times g$ subblock of $\mathbf{H}$, then

[TABLE]

where $\mathbf{N}$ is the nilpotent matrix with ones on its superdiagonal and zeros everywhere else. Then $4\pi\mathrm{i}\tilde{\mathbf{A}}^{-1}\mathbf{H}_{g}\mathbf{L}^{-1}$ is a Riemann matrix (e.g., it has a negative definite real part). Since this matrix is real, we conclude that $\mathbf{H}_{g}$ is invertible.

Furthermore,

[TABLE]

This follows from the fact that this summation can be written as

[TABLE]

where $\gamma$ is a clockwise oriented curve encircling $[a_{1},b_{g+1}]$. This integral can be evaluated by taking the residue at infinity using $R(z)\sim z^{g+1}$ as $z\to\infty$. Thus, we have that

[TABLE]

where $\mathbf{0}_{g}$ is the zero vector of length $g$, $\mathbf{1}_{g}$ is a vector of length $g$ of all ones, and $\mathbf{H}_{g+1,1:g}$ is the first $g$ columns of the $(g+1)$th row of $\mathbf{H}$. This new matrix has invertible diagonal blocks, so it must be invertible, and thus $\mathbf{H}$ is also invertible. ∎

By construction, we have that

[TABLE]

Now, consider the transformation

[TABLE]

Noting that if $\mathbf{S}_{n}^{+}(z)=\mathbf{S}_{n}^{-}(z)\mathbf{K}(z)$ on a contour, then

[TABLE]

we find that for $z\in(a_{j+1},b_{j})$,

[TABLE]

meaning that the jump condition on the gaps has been eliminated. Applying our transformation to the remaining jump conditions, we use the fact that the asymptotic conditions at the endpoints are unaffected since $|\mathfrak{h}_{n}|$ is bounded by [1, Lemma 7.2.2]. We arrive at the new Riemann–Hilbert problem

[TABLE]

that we will solve numerically. As shorthand, we define

[TABLE]

We include example plots of the curves and jump conditions in this problem in Figure 3.

3.3. Recurrence coefficients

Recurrence coefficients for the orthogonal polynomials corresponding to the weight function (1) are directly related to the first-order asymptotic behavior (i.e., the residue at infinity) of the unique solution (5) to the original problem (6). If

[TABLE]

as $z\to\infty$, then the recurrence coefficients (2) are given by

[TABLE]

for $n=0,1,2,\ldots$ [9]. We now use this to derive equivalent formulae based on the first-order asymptotic behavior of (14) by working backwards through the transformations. Define $\mathfrak{h}_{n}^{(1)}$ by

[TABLE]

Then,

[TABLE]

From this, if

[TABLE]

as $z\to\infty$, then

[TABLE]

implying $\mathbf{S}_{n}^{(1)}=\tilde{\mathbf{S}}_{n}^{(1)}-\mathfrak{h}_{n}^{(1)}\sigma_{3}$.

Now, we obtain the first-order asymptotics of $\mathbf{Z}_{n}$ (7) from $\mathbf{S}_{n}$ (11). Recall that

[TABLE]

Then, the asymptotic expansion of $\mathfrak{g}$ (10) gives

[TABLE]

Expanding

[TABLE]

gives

[TABLE]

and therefore

[TABLE]

Finally, the lensing of the problem in Step 1 does not affect the asymptotic behavior at infinity, so

[TABLE]

yielding new recurrence formulae

[TABLE]

that we will use to recover the recurrence coefficients when numerically solving the problem (14).

As $n\to\infty$, the jump conditions on $C_{j}$ tend exponentially to the identity matrix, so the FIK Riemann–Hilbert problem can be considered with just the jump conditions on the intervals, asymptotically. While formulae for the recurrence coefficients in this asymptotic case are known [12, Appendix B], their evaluation is made difficult due to the presence of theta functions composed with an Abel map. Thus, a numerical method for solving this form of the FIK Riemann–Hilbert problem not only allows for efficient evaluation of these rather complicated formulae, but can potentially be utilized to compute theta functions. Related to this, one could choose the functions $h_{j}$ so that $w_{j}(z)=w_{k}(z)$ for all $j,k$. In this case, individual circles around each interval $[a_{j},b_{j}]$ could be replaced with a large circle around $[a_{1},b_{g+1}]$. We will explore this further in future work.

4. Numerical methods

4.1. Numerically computing Cauchy integrals

We compute orthogonal polynomial coefficients via (18) by solving the Riemann–Hilbert problem (14) numerically. Doing so will heavily rely on the computation of Cauchy integrals and, specifically, on computing Cauchy integrals of Chebyshev polynomials on an interval. Recalling that the weight function of Chebyshev-$\mathcal{T}$ on $\mathbb{I}=[-1,1]$ is given by $w_{\mathcal{T}}(x)=\frac{1}{\pi\sqrt{x+1}\sqrt{1-x}}$ and that $\mathcal{T}_{0}=1$,

[TABLE]

Similarly,

[TABLE]

Theorem 4.1**.**

Let $C^{\square}_{k}(z)=\mathcal{C}_{\mathbb{I}}[\square_{k}w_{\square}](z)$ where $\square$ denotes the kind of Chebyshev polynomial $\mathcal{T},U,V,W$. Then, for $k\geq 1$,

[TABLE]

where $J^{-1}_{+}$ is a right inverse of the Joukowsky map, mapping $\mathbb{C}\setminus[-1,1]$ to the interior of the unit circle, given by

[TABLE]

Proof.

We proceed by induction by first considering the base case $k=1$ separately for each kind of Chebyshev polynomial. Using the recurrence relation (3) that $C^{\square}_{k}$ must satisfy, we can directly compute the following.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Now,

[TABLE]

Using this for the second base case $k=2$,

[TABLE]

As an inductive step, we let $k\geq 3$ and find that our recurrence gives

[TABLE]

completing the proof for Chebyshev-$\mathcal{T}$. For $U,V$, and $W$, see that the theorem is trivially true for $k=0$, so we can use $k=0,1$ as our base cases. Letting $\square$ denote $U,V$, or $W$, we let $k\geq 2$ and find from the recurrence that

[TABLE]

completing the proof for all cases. ∎

These integrals can be scaled and shifted appropriately to compute $\mathcal{C}_{[a,b]}[fw]$ for an arbitrary interval $[a,b]$. If we let $\tilde{\mathcal{T}}_{k}(x;[a,b])=\mathcal{T}_{k}(M^{-1}(x;[a,b]))$ with weight function $\tilde{w}_{\mathcal{T}}(x;[a,b])=\frac{1}{\pi\sqrt{x-a}\sqrt{b-x}}$, $\tilde{U}_{k}(x;[a,b])=U_{k}(M^{-1}(x;[a,b]))$ with weight function $\tilde{w}_{U}(x;[a,b])=\frac{2}{\pi}\left(\frac{2}{b-a}\right)^{2}\sqrt{x-a}\sqrt{b-x}$, $\tilde{V}_{k}(x;[a,b])=V_{k}((M^{-1}(x;[a,b])))$ with weight function $\tilde{w}_{V}(x;[a,b])=\frac{2}{\pi(b-a)}\frac{\sqrt{x-a}}{\sqrt{b-x}}$, and $\tilde{W}_{k}(x;[a,b])=W_{k}(M^{-1}(x;[a,b]))$ with weight function $\tilde{w}_{W}(x;[a,b])=\frac{2}{\pi(b-a)}\frac{\sqrt{b-x}}{\sqrt{x-a}}$ where

[TABLE]

then

[TABLE]

It is important to determine the leading-order asymptotic behavior of these integrals for our later calculations. We have

[TABLE]

where $\square$ denotes any of the four kinds of Chebyshev polynomial as $z\to\infty$. This follows from the above normalizations.

4.2. Numerically computing $\mathfrak{g}$

To determine $\mathfrak{g}$, we first obtain the coefficients of $Q_{g}(z)$ via solving the linear system (8) by numerically computing the relevant integrals via a Chebyshev-$\mathcal{T}$-based quadrature rule222We can easily do this by getting the first Chebyshev coefficient of $z^{j}\sqrt{z-a_{k}}\sqrt{b_{k}-x}/R(z)$ in a mapped Chebyshev-$\mathcal{T}$ series via the discrete cosine transform..

Lemma 4.2** (Plemelj Lemma, [29], Section 17).**

Let $\Gamma\subset\mathbb{C}$ and $f$ be Hölder continuous on $\Gamma$, then for $z\in\Gamma$ that is not an endpoint of $\Gamma$,

[TABLE]

Motivated by [39, Section 5.5], we use Lemma 4.2 to rewrite $\mathfrak{g}^{\prime}$ as a Cauchy integral. We have that $R^{+}(z)=-R^{-}(z)$ for $z\in(a_{j},b_{j})$, and $Q_{g}$ is of lower degree than $R$, so $\mathfrak{g}^{\prime}$ solves the scalar Riemann–Hilbert problem

[TABLE]

Thus, by Lemma 4.2,

[TABLE]

Of course, this implies that

[TABLE]

where the integration path is taken to be a straight line between $a_{1}$ and $z$. We approximate $\frac{Q_{g}(z)}{R^{+}(z)}\sqrt{z-a_{j}}\sqrt{z-b_{j}}$ to high-precision as a Chebyshev-$\mathcal{T}$ series on $[a_{j},b_{j}]$ using the discrete cosine transform. Letting $w_{\mathcal{T},j}(x)=\tilde{w}_{\mathcal{T}}(x;[a_{j},b_{j}])$, $\mathcal{T}_{j,k}(x)=\tilde{\mathcal{T}}_{k}(x;[a_{j},b_{j}])$ and dropping the lower bound of integration that will be accounted for later, it thus suffices to consider computing the anti-derivatives

[TABLE]

where $\alpha_{j,k}$ is the $k$th Chebyshev-$\tilde{\mathcal{T}}$ coefficient of $\frac{Q_{g}(z)}{R^{+}(z)}\sqrt{z-a_{j}}\sqrt{z-b_{j}}$ on $[a_{j},b_{j}]$. For $k>0$, we swap the order of integration

[TABLE]

Let $f(s)=\mathcal{T}_{j,k}(s)w_{\mathcal{T},j}(s)$. Using integration by parts,

[TABLE]

where

[TABLE]

By orthogonality,

[TABLE]

for $k>0$, and333We use $\diamond$ notation to specify a function without naming it. For example, we use $1/\diamond$ to denote the function $f(x)=1/x.$

[TABLE]

Now, [13, Eq. 18.9.22] gives that

[TABLE]

where $\mathcal{T}_{j}(x)=\sqrt{2}T_{j}(x)$ in our convention. Substituting this in and integrating both sides gives

[TABLE]

and

[TABLE]

Using the change of variables $t=M_{j}^{-1}(s)$ where $M_{j}^{-1}(s)=M^{-1}(s;[a_{j},b_{j}])$,

[TABLE]

Recall that by definition

[TABLE]

and

[TABLE]

so that

[TABLE]

where we use the same notation $w_{U,j}(x)=w_{U}(x;[a_{j},b_{j}])$, $U_{j,k}(x)=\tilde{U}_{k}(x;[a_{j},b_{j}])$. Thus,

[TABLE]

for $k>0$. For $k=0$,

[TABLE]

which we have computed using the exact formulae for the Cauchy integrals of Chebyshev-$\mathcal{T}$ polynomials with their respective weight functions. With a change of variables $t=M(z^{\prime})$ as before, we find [31]

[TABLE]

Set

[TABLE]

where $m_{j}$ is the number of terms we use in each Chebyshev series. Then, $\mathfrak{g}$ is approximated by

[TABLE]

A numerical instability can arise when computing $\Phi(a_{1})$ with this approach due to the non-Lipschitz nature of $J^{-1}_{+}(z)$ at $z=-1$. This can be circumvented by replacing $J^{-1}_{+}(-1)=-1$ at each instance in the evaluation of $\Phi(a_{1})$ that it occurs.

4.2.1. Asymptotic behavior

We compute the logarithmic444This does not appear in our recurrence formulae but is necessary in (31) below. Additionally, $|\mathfrak{c}|$ is commonly referred to as the capacity of $\cup_{j}[a_{j},b_{j}]$ [33]. and first order terms $\mathfrak{c}$ and $\mathfrak{g}_{1}$ in the asymptotic expansion of $\mathfrak{g}$ (10). Recall that $\mathfrak{g}$ can be numerically computed as the Cauchy integrals of Chebyshev-$U$ polynomials and the logarithm of the inverse Joukowsky map, so it suffices to compute the asymptotics of these functions.

The only term in our expansion (21) that contains a logarithmic term is (20). We have

[TABLE]

as $z\to\infty$, so that

[TABLE]

The only other contribution to the $\operatorname{O}(z^{-1})$ term comes from

[TABLE]

by (19). This gives

[TABLE]

The first order term can be approximated as

[TABLE]

4.3. Numerically computing $\mathfrak{h}_{n}$

To numerically compute $\mathfrak{h}_{n}$ we first compute $\Delta_{1},\ldots,\Delta_{n}$ according to (9) by numerically evaluating $\mathfrak{g}^{+}(z)+\mathfrak{g}^{-}(z)$ on the gaps $(b_{j},a_{j+1})$. We then obtain $A_{1}(n),\ldots,A_{g+1}(n)$ in the same manner as we obtained the coefficients for $\mathfrak{g}$. We do this by numerically computing the integrals via a Chebyshev-$\mathcal{T}$-based quadrature rule and solving the system (13).

Once the coefficients are known, we can evaluate $\mathfrak{h}_{n}$ in the complex plane in accordance with (12) by developing a method to compute the relevant Cauchy integrals. We do so by utilizing the formulae of Theorem 4.1. Since we know $R(z)$ and our coefficients, we need only numerically evaluate

[TABLE]

for $j=1,\ldots,g+1$, $\ell=1,\ldots,g$. To do this, let $a=b_{j}$, $b=a_{j+1}$ or $a=a_{j}$, $b=b_{j}$ for some $j$. If we consider a truncated Chebyshev series

[TABLE]

Then,

[TABLE]

Of course, these Cauchy integrals can be computed via our exact formulae in Theorem 4.1, so we can approximate $\mathfrak{h}_{n}$ by computing a truncated series and taking linear combinations of Cauchy integrals as appropriate.

4.3.1. First-order asymptotic behavior

Now, we compute the first order term $\mathfrak{h}_{n}^{(1)}$ in (17). This is the $\operatorname{O}(z^{-g-2})$ term of the Cauchy integrals in (12) since $R(z)\sim z^{g+1}$ as $z\to\infty$. To compute this, we use the asymptotics of Cauchy integrals to find

[TABLE]

To compute this term numerically, we use a Chebyshev-$\mathcal{T}$-based quadrature rule on each interval to evaluate the integrals.

4.4. Numerically solving the Riemann–Hilbert problem and recovering coefficients

Now that we can compute $\mathfrak{g}$ and $\mathfrak{h}_{n}$ reliably in the complex plane, we develop a method to solve the Riemann–Hilbert problem (14) numerically. We do so by developing methods to solve smaller and simpler problems that we then combine into a block linear system.

4.4.1. Numerically solving scalar Riemann–Hilbert problems on a circle

As a first step, we numerically solve a Riemann–Hilbert problem with a jump condition occurring only on the unit circle. Letting $\mathbb{U}$ denote the unit circle in the complex plane oriented counterclockwise, consider the problem of finding a function $\Phi:\mathbb{C}\setminus\mathbb{U}\to\mathbb{C}$ that satisfies

[TABLE]

To find a numerical solution, we make the ansatz

[TABLE]

for some unknown function $u$. This ansatz can be rigorously justified, see [39, Section 2.5]. Plugging this into the jump condition in (22),

[TABLE]

Now, we let $u$ be given by a Laurent series

[TABLE]

A simple residue calculation yields that

[TABLE]

To form a numerical solution, we truncate our Laurent series so that its terms range from $-N$ to $N$ and consider $2N+1$ collocation points $z_{0},\ldots,z_{2N}$ evenly distributed on $\mathbb{U}$ 555For instance, given by $z_{k}=\mathrm{e}^{2\pi ik/(2N+1)}$.. We then determine coefficients $c_{-N},\ldots,c_{N}$ such that (23) is satisfied at all collocation points:

[TABLE]

for $k=0,\ldots,2N$. To solve for the coefficients, we write (25) as a linear system and solve for $c_{-N},\ldots,c_{N}$. We then use (24) to evaluate our numerical solution $\Phi(z)$ at any $z\in\mathbb{C}\setminus\mathbb{U}$. One can of course generalize this to numerically solve a scalar Riemann–Hilbert problem on any circle by scaling and shifting appropriately [30].

4.4.2. Numerically solving scalar Riemann–Hilbert problems on an interval

Now, we numerically solve the Riemann–Hilbert problem of finding a function $\Phi:\mathbb{C}\setminus\mathbb{I}\to\mathbb{C}$ that satisfies

[TABLE]

where $\mathbb{I}=[-1,1]$. We again make the ansatz

[TABLE]

but need a different approach to determine $u$. As with most numerical analysis problems on the unit interval, we cannot safely use a monomial basis and equally-spaced points as in the problem on the unit circle. Instead, we need to use an orthogonal polynomial basis, meaning that we will need to use the Cauchy integrals of orthogonal polynomials when enforcing our jump condition of the form (23) but on $\mathbb{I}$. We choose $2N+1$ collocation points to be more concentrated near the endpoints as the roots of Chebyshev-$\mathcal{T}$, i.e., we set

[TABLE]

If we use an orthogonal polynomial basis $p_{0},\ldots,p_{2N}$ with associated weight function $w$, then we find coefficients $d_{0},\ldots,d_{2N}$ such that the jump condition holds at all collocation points:

[TABLE]

for $k=0,\ldots,2N$. This is a linear system and can (potentially) be solved for $d_{0},\ldots,d_{2N}$ by computing $\mathcal{C}_{\mathbb{I}}^{+}[p_{j}w](z_{k}),\mathcal{C}_{\mathbb{I}}^{-}[p_{j}w](z_{k})$ as described in accordance with Theorem 4.1 if we choose to use a Chebyshev basis for $p_{j}$, $j=0,\ldots,2N$. Of course, this problem can be shifted and scaled to an arbitrary interval $[a,b]$ [30].

While solving this linear system appears straightforward, a potential difficulty arises in choosing the orthogonal polynomial basis used to build the numerical solution. If the correct basis is not chosen, the method will not converge. The Cauchy integrals of different orthogonal polynomial bases have different asymptotic behavior near the endpoints $z\to\pm 1$, and therefore different choices of bases will produce functions that behave very differently. If one wishes to pose a Riemann–Hilbert problem with a unique solution, additional asymptotic conditions must be prescribed at the endpoints. Then, the correct orthogonal polynomial basis associated with that asymptotic condition must be selected.

4.4.3. Combining methods and solving the orthogonal polynomial problem

Consider solving a Riemann–Hilbert problem on the contour $\Gamma=\Gamma_{1}\cup\Gamma_{2}$ where $\Gamma_{1},\Gamma_{2}$ are disjoint contours: Find a function $\Phi:\mathbb{C}\setminus\Gamma\to\mathbb{C}$ that satisfies

[TABLE]

Define the projection operator $P_{j}$ by

[TABLE]

The ansatz becomes

[TABLE]

where $u_{j}=P_{j}u$ for some unknown $u$. Noting that $\mathcal{C}_{\Gamma}=\mathcal{C}_{\Gamma_{1}}+\mathcal{C}_{\Gamma_{2}}$, the jump condition becomes

[TABLE]

for $j,k\in\{1,2\}$, $j\neq k$. This gives the block system

[TABLE]

where $P_{j}C_{\Gamma_{k}}^{+}=P_{j}C_{\Gamma_{k}}^{-}=P_{j}C_{\Gamma_{k}}$ when $j\neq k$. We expand $u_{1},u_{2}$ in the usual way666It is expanded using either a Laurent series or an appropriate orthogonal polynomial basis depending on if $\Gamma_{1},\Gamma_{2}$ are intervals or circles. and choose collocation points for each on $\Gamma_{1}$ and $\Gamma_{2}$, respectively. Then, we enforce the jump conditions at the collocation points as a block linear system via the matrix structure (27) and solve for the basis coefficients for both $u_{1}$ and $u_{2}$. A block system can be formulated analogously when $\Gamma$ is composed of an arbitrary number of disjoint contours. In this case, the $(j,k)$th entry in the block matrix is given by $P_{j}C_{\Gamma_{k}}^{+}-f_{j}P_{j}C_{\Gamma_{k}}^{-}$.

We can extend these methods to a matrix Riemann–Hilbert problem for $\mathbf{\Phi}:\mathbb{C}\setminus\Gamma\to\mathbb{C}^{2\times 2}$

[TABLE]

If one considers the ansatz

[TABLE]

where $\mathbf{U}$ is now matrix-valued and the Cauchy operator is taken element-wise. We find a system of the form

[TABLE]

We solve this using our methods for scalar problems by forming a block system. Solving a matrix problem on multiple domains will require forming a block system of block systems [30]. The equations for $\mathbf{U}_{11},\mathbf{U}_{12}$ and $\mathbf{U}_{21},\mathbf{U}_{22}$ decouple, so $\mathbf{U}$ can be solved for row-wise.

The pieces are now in place to numerically solve the Riemann–Hilbert problem (14). As a simplifying assumption, we take each $C_{j}$ in Figure 2 to be a circle centered at $(a_{j}+b_{j})/2$ with a diameter greater than $b_{j}-a_{j}$777Heuristically, we find that the diameter should be at least $5(b_{j}-a_{j})/4$ to avoid stability issues. but small enough so that $C_{j}$ and $C_{k}$ do not intersect when $j\neq k$888To be fully robust, one should take $C_{j}$ to be a Bernstein ellipse to which $1/h_{j}(z)$ is analytically continuable, but this requires a method to numerically solve Riemann–Hilbert problems on such an ellipse.. Then, (14) can be solved via a block linear system according to our methods for circles, intervals, multiple domains, and matrix problems by accounting for shifting and scaling things to the unit interval and circle as appropriate. The only remaining difficulty lies in choosing orthogonal polynomial bases so that the asymptotic conditions at the endpoints are satisfied.

On each interval $[a_{j},b_{j}]$, we choose orthogonal polynomial bases based on the endpoint asymptotic conditions from (14):

[TABLE]

Because the weight functions have asymptotic behavior

[TABLE]

we choose bases as follows:

• •

If $\alpha_{j},\beta_{j}=-1$, then we use $\mathcal{C}_{[a_{j},b_{j}]}[U_{j,k}w_{U,j}]$ in the first column of the matrix ansatz $\mathbf{U}$ (28) and $\mathcal{C}_{[a_{j},b_{j}]}[\mathcal{T}_{j,k}w_{\mathcal{T},j}]$ in the second column of the ansatz.

• •

If $\alpha_{j},\beta_{j}=1$, then we use $\mathcal{C}_{[a_{j},b_{j}]}[\mathcal{T}_{j,k}w_{\mathcal{T},j}]$ in the first column and $\mathcal{C}_{[a_{j},b_{j}]}[U_{j,k}w_{U,j}]$ in the second column.

• •

If $\alpha_{j}=1,\beta_{j}=-1$, then we use $\mathcal{C}_{[a_{j},b_{j}]}[W_{j,k}w_{W,j}]$ in the first column and $\mathcal{C}_{[a_{j},b_{j}]}[V_{j,k}w_{V,j}]$ in the second column.

• •

If $\alpha_{j}=-1,\beta_{j}=1$, then we use $\mathcal{C}_{[a_{j},b_{j}]}[V_{j,k}w_{V,j}]$ in the first column and $\mathcal{C}_{[a_{j},b_{j}]}[W_{j,k}w_{W,j}]$ in the second column.

4.4.4. Recovering coefficients

Once a numerical solution to (14) is obtained, we compute $\tilde{\mathbf{S}}_{n}^{(1)}$ by summing the first order terms of each approximate function in the block system. More explicitly, if the $(\ell,m)$ entry of $\tilde{\mathbf{S}}_{n}$ has a Laurent series representation

[TABLE]

on each $C_{j}$ where $q_{j}=(a_{j}+b_{j})/2$ and $r_{j}$ denotes the radius of $C_{j}$ and an orthogonal polynomial series representation

[TABLE]

on each $[a_{j},b_{j}]$, then

[TABLE]

where the scaling follows from (19) if $p_{k}$ and $w$ are a shifted and scaled Chebyshev polynomial and its weight, respectively. This combined with our computed first order contributions of $\mathfrak{g}$ and $\mathfrak{h}_{n}$ from earlier allows us to recover our recurrence coefficients from (18) as

[TABLE]

5. Examples and applications

5.1. The single interval case

As a first example, we consider our methods applied to a single interval for which recurrence coefficients are known asymptotically. For simplicity, we consider this interval to be $[-1,1]$. It turns out that the functions $\mathfrak{g}$ and $\mathfrak{h}_{n}$ can be computed exactly in this case, so solving the Riemann–Hilbert problem (14) numerically becomes simpler. First, we observe that there are no gaps between intervals on which jump conditions can be induced, so the function $\mathfrak{h}_{n}$ is identically zero for all $n$. Furthermore, determining $\mathfrak{g}^{\prime}$ becomes trivial, and the exponential of $\mathfrak{g}$ can be computed explicitly as [24]

[TABLE]

and the weight function is given by

[TABLE]

and $h_{1}$ satisfies the usual assumptions of being positive on $[-1,1]$ and analytic in a neighborhood of it [25]. We additionally find that the first-order correction $\mathfrak{g}_{1}=0$ in this case, meaning that all correction terms in the recurrence coefficient formulae (18) drop out.

As with the multiple interval problem, we choose $C_{1}=\{z:|z|=5/4\}$, restricting our choice of $h$ to functions that are positive on $[-1,1]$ and whose reciprocals are analytic in an open set containing $\bar{D}_{1}$. For our example, we choose

[TABLE]

and $\alpha_{1}=\beta_{1}=1$ to obtain a Chebyshev-$U$-like weight. We first compute the first 51 coefficients of these orthogonal polynomials using 160 collocation points on $C_{1}$ and 16 collocation points on $[-1,1]$. We should expect these coefficients $a_{j},b_{j}$ to tend to the Chebyshev-$U$ coefficients ($a_{j}=0,b_{j}=1/2$) exponentially as $j\to\infty$ [16]. In Figure 4, we plot these differences $|a_{j}|$ and $|b_{j}-1/2|$ for the first 51 computed coefficients that exhibit the expected behavior.

Now, we present an alternative Householder method for computing recurrence coefficients, to which the accuracy of our method is compared. Consider a discrete $N$-point approximation to the weight $w$ with nodes $x_{1},\ldots,x_{N}$ and weights $w_{1},\ldots,w_{N}$, i.e.,

[TABLE]

Define

[TABLE]

Then, from [15, Eq. 2.2.24], if the nodes and weights are selected according to the Gauss quadrature rule for $w$, we have a factorization of the form

[TABLE]

where $\mathbf{e}_{1}$ is the first unit vector and $\mathbf{J}_{N}(w)$ is the $N\times N$ principal subblock of $\mathbf{J}(w)$. If one can compute $\mathbf{Q}$, this allows for the recovery of $\mathbf{J}_{N}(w)$; the Householder algorithm [22] is a stable way to do this. Due to the complexity of the Householder algorithm, this approach requires $\operatorname{O}(N^{3})$ arithmetic operations to compute the first $N$ recurrence coefficients. An algorithm based on Givens rotations with $\operatorname{O}(N^{2})$ arithmetic operations is presented in [19]. An important distinction between these two methods and the Stieltjes procedure is that these two methods require a priori knowledge of the desired recurrence length $N$, and altering this a posteriori requires restarting the entire algorithm.

In Figure 5, we plot the error at each of the first 51 coefficients for various numbers of collocation points on the $[-1,1]$, always choosing the number of collocation points on $C_{1}$ to be 10 times the number on $[-1,1]$.

The method is able to compute the recurrence coefficients $a_{j},b_{j}$ starting at any choice of $j$ without first computing $a_{k},b_{k}$ for $k<j$. To demonstrate this, we include a similar plot showing the errors in computing $a_{j}$ and $b_{j}$ for $j=1000,\ldots,1020$ for various numbers of collocation points in Figure 6.

5.2. Multiple intervals

Now, the accuracy of our method is demonstrated in computing with orthogonal polynomials on multiple intervals. In all examples, we will take $h_{j}(x)=1$ for all $j$, for simplicity, and compare our computed coefficients against those computed by the Householder method described in Section 5.1.

In Figure 7, we consider orthogonal polynomials defined on $[-1.8,-1]\cup[2,3]$ with a Chebyshev-$\mathcal{T}$-like weight on each ($\alpha_{j}=\beta_{j}=-1$ for $j=1,2$).

Furthermore, in Figure 8, we include a timing comparison of our method to an optimized $\operatorname{O}(N^{2})$ algorithm given as RKPW in [19] and lanczos.m in [15].

Since we are concerned with asymptotic timings for this comparison, our method is applied with only the jump conditions on the intervals as mentioned in Section 3.3. As is also mentioned in this section, there exist classes of weight functions for which the circles $C_{j}$ are not needed for any polynomial degree, and we will explore this further in future work. We observe that our method begins to outperform the RKPW method once around 3200 pairs of recurrence coefficients need to be computed; however, since the RKPW method is heavily optimized while our method is coded suboptimally, we expect the threshold to be much lower with a better implementation. The RKPW method also has the aforementioned weakness in that it requires a priori knowledge of the desired recurrence length $N$, and altering this a posteriori requires potentially restarting the entire algorithm. Thus, in practice, one should consider using the RKPW method to compute the first few coefficients and our method asymptotically.

In Figure 9, we consider orthogonal polynomials defined on $[0.1,1.1]\cup[2,3]\cup[3.5,4]$ with a Chebyshev-$V$-like weight on each ($\alpha_{j}=1,\beta_{j}=-1$ for $j=1,2,3$). Finally, in Figure 10, we consider orthogonal polynomials defined on $[-3.2,-2.2]\cup[0.1,1.1]\cup[2,3]\cup[3.5,4]$ with a Chebyshev-$\mathcal{T}$-like weight on $[-3.2,-2.2]$, a Chebyshev-$U$-like weight on $[0.1,1.1]$, a Chebyshev-$V$-like weight on $[2,3]$, and a Chebyshev-$W$-like weight on $[3.5,4]$ ($\alpha_{1}=\beta_{1}=-1$, $\alpha_{2}=\beta_{2}=1$, $\alpha_{3}=1,\beta_{3}=-1$, $\alpha_{4}=-1,\beta_{4}=1$). In each example, we use the same number (2, 4, 8, and 16) of collocation points on each interval and 10 times that number on each circle $C_{j}$ and compute the first 51 recurrence coefficients. In practice, one should consider using a different number of collocation points on each interval and varying the ratio of collocation points on circles to points on intervals to optimize runtime.

Heuristically, our method as described encounters stability issues once $g$ reaches 4 or 5; however, 6 digits of accuracy is maintained until about $g=11$. This is because the linear systems (8) and (13) become ill-conditioned as $g$ increases. If one wishes to consider larger $g$, these instabilities can be circumvented by instead building linear systems with the basis in [4, Eq. 28].

In all figures, we observe rapid convergence to the “ground truth” coefficients $a_{j},b_{j}$ for any value of $j$. Errors on the order of machine precision are always achieved when we use 16 collocation points on each interval and 160 points on each circle. Oftentimes, less collocation points are needed to achieve this accuracy.

5.3. The Toda lattice

The semi-infinite Toda lattice as considered by Deift, Li, and Tomei [10] (see also [28]) with boundary condition $b_{-1}(t)=0$ is defined by semi-infinite operators

[TABLE]

using the relation

[TABLE]

where $[\cdot,\cdot]$ denotes the standard matrix commutator. Given an initial condition $\mathbf{X}(0)=\mathbf{J}(w)$ corresponding to the Jacobi operator to some orthogonal polynomial weight $w$, the time-dependent $a_{j}(t),b_{j}(t)$ are given by the respective recurrence coefficients $a_{j},b_{j}$ of the orthogonal polynomials with weight $\mathrm{e}^{tx}w(x)$ [10] (see also [35]), i.e., if $\mathbf{X}(0)=\mathbf{J}(w)$, then $\mathbf{X}(t)=\mathbf{J}(w\mathrm{e}^{t\diamond})$.

Since our method allows for weights to be scaled by positive analytic functions, we are able to numerically evolve the Toda lattice for initial conditions corresponding to the various orthogonal polynomial weights in our arsenal; however, the scaling of $\mathrm{e}^{tx}$ grows exponentially in the complex plane, so we are only able to evolve the Toda lattice to a finite time horizon before encountering precision issues.

As a first example, we consider an initial condition of

[TABLE]

corresponding to a Chebyshev-$U$ weight on $[-1,1]$. Choosing to use 120 collocation points on $C_{1}$ and 20 on $[-1,1]$, we compute the first 11 recurrence coefficients with $h_{1}(x)=\mathrm{e}^{tx}$ for times ranging from $t=0$ to $t=11.5$ and plot the results in Figure 11. Numerical blowup is observed shortly after time $t=11.5$.

Now, we compute the time-evolution of the Toda lattice with initial condition $\mathbf{X}(0)$ corresponding to the recurrence coefficients of orthogonal polynomials with a Chebyshev-$\mathcal{T}$-like weight on $[-3,-2]\cup[2,3]$. Numerically, we observe that $a_{j}(0)=0$ for all $j$ and as $j\to\infty$, $b_{j}(0)$ alternates between $0.5$ and $2.5$ as displayed in Figure 12.

We evolve the Toda lattice numerically for times ranging from $t=0$ to $t=4.5$ in the same manner as before999By taking $h_{j}(x)=\mathrm{e}^{tx}$ for each $j$ and computing the first 11 recurrence coefficients at each time with 120 collocation points on each $C_{j}$ and 20 on each interval, we approximate the evolution of the Toda lattice. and plot the results in Figure 13.

Additionally, we emphasize that our method allows us to compute any recurrence coefficient without first computing the others, meaning that the evolution of $a_{j}(t),b_{j}(t)$ can be computed for any $j$. We demonstrate this by plotting $a_{116}(t),b_{116}(t)$ in Figure 14.

Finally, we consider the time-evolution of the Toda lattice with initial condition $\mathbf{X}(0)$ corresponding to the recurrence coefficients of orthogonal polynomials with a Chebyshev-$W$-like weight on $[0.1,1.1]\cup[2,3]\cup[3.5,4]$. Now, the initial condition as displayed as in Figure 15 appears almost noisy.

We evolve the Toda lattice numerically for times ranging from $t=0$ to $t=5.7$ in the same manner as before and plot the results in Figure 16.

5.4. Building an approximation to $1/x$

As a final example, we demonstrate that our procedure can be used to approximate functions on disconnected domains. Specifically, we consider an approximation to the function $1/x$. If the orthonormal polynomials corresponding to the normalization of the weight (1) are given by $p_{0},p_{1},\ldots$, then a $p_{j}$-series expansion for $1/x$, for $x\in\Sigma$, is given by

[TABLE]

where $\Sigma=\bigcup_{j=1}^{g+1}[a_{j},b_{j}]$, and $\eta=\int_{\Sigma}w(x)\mathrm{d}x$. Allowing for unnormalized weight functions mandates the inclusion of the $\frac{1}{\eta}$ constant that would not be present if $w$ was a probability measure or if $p_{0},p_{1},\ldots$ corresponded to the unnormalized weight. Our definition of $p_{0},p_{1},\ldots$ allows for the Cauchy integrals of these polynomials to be computed without also computing $\eta$. Thus, we need only to evaluate the Cauchy integral of our polynomials at zero to build an approximation. To do this, we recall that the (1,2) entry of the true solution (5) to the original Riemann–Hilbert problem (6) is the Cauchy integral of the monic orthogonal polynomial, meaning that

[TABLE]

where $\gamma_{n}$ is the coefficient on the $x^{n}$ term in $p_{n}(x)$. [9] gives that

[TABLE]

with $\gamma_{0}=1$, so this term can be easily computed if we are already computing the recurrence coefficients. To obtain this from a computed solution $\tilde{\mathbf{S}}_{n}(z)$, we undo our series of transformations. If we assume that $z\notin\Sigma$ and furthermore that $z\notin\bar{D}_{j}$ where $D_{j}$ is the interior of $C_{j}$ for any $j$, then we need not worry about undoing the lensing step. In fact, the second column of $\mathbf{Y}_{n}$ is unaffected by the lensing step, so this need not be accounted for when computing $\mathcal{C}_{\Sigma}\left[p_{n}w\right](z)$, regardless of whether $z\in\bar{D}_{j}$. Since the remaining steps are diagonal, we can compute

[TABLE]

For our purposes, we only need to evaluate this at $z=0$, so this can easily be done in parallel with computing the recurrence coefficients by computing $\mathfrak{g}(0)$ once at the beginning and computing $\mathfrak{h}_{n}(0)$ each time we solve the Riemann–Hilbert problem (14).

Using (30) and (31), we build an approximation to $1/x$ on $\Sigma=[-4,-3]\cup[-2,-1]\cup[2,3]$ as well as $\Sigma=[1,2.3]\cup[3,4]\cup[4.5,6.1]$ where each Riemann–Hilbert problem is numerically solved with 120 collocation points on each $C_{j}$ and 20 points on each interval. To check accuracy, we consider an evenly spaced grid by 0.01 $x_{1},x_{2},\ldots$ on each domain $\Sigma$ and consider the max-norm error

[TABLE]

where $N$ is the number of polynomials that we use. In Figure 17, we plot this computed error for various numbers of terms in the expansion (30).

It is also possible for us to determine the rate of convergence for our orthogonal polynomial approximation. From [12], we have that if [math] is outside the deformation region of our Riemann–Hilbert problem,

[TABLE]

where $\delta_{n}$ depends on $\prod_{j=0}^{n-1}\frac{1}{b_{j}\mathfrak{c}}$,and other constants determined by the asymptotic behavior of $\mathbf{Y}_{n}$ at infinity [12]. By [12, Equations B.3, B.12],

[TABLE]

as $n\to\infty$. This implies that $\delta_{n}$ is bounded from above, so we have that

[TABLE]

where $c$ and $c_{1}$ are constants. This implies that our approximation converges uniformly, at worst, at a geometric rate given by $|\mathrm{e}^{-\mathfrak{g}(0)}|$. We believe $\delta_{n}$ to also be bounded below which would imply that convergence is truly given by this rate. This means that in practice, one can compute $\mathfrak{g}(0)$ and the error at an initial approximation that will then allow them to estimate the number of terms necessary to reach a desired accuracy.

6. Conclusion and future work

This work presents a novel numerical method for computing orthogonal polynomials that are orthogonal on multiple, disjoint intervals, where analytical formulae are currently unavailable. By leveraging the Fokas–Its–-Kitaev Riemann–-Hilbert representation of these polynomials, our approach achieves $\operatorname{O}(N)$ complexity for computing the first $N$ recurrence coefficients. Moreover, our method enables efficient pointwise evaluation of the polynomials and their Cauchy transforms across the entire complex plane.

The efficiency of our method is currently limited for low degree polynomials by the jump conditions on the circles $C_{j}$ around each interval; however, there exist functions $h_{j}$ such that the analytic continuations of $w_{j},w_{k}$ satisfy $w_{j}(z)=w_{k}(z)$ for all $j,k$ and individual circles around each interval $[a_{j},b_{j}]$ could be replaced with a large circle around $[a_{1},b_{g+1}]$, and then removed entirely. In such a regime, our method is competitive (and better asymptotically) when compared with the $\operatorname{O}(N^{2})$ method presented in [19, 15]. In future work, we will explore these classes of weight functions and apply their respective orthogonal polynomials to indefinite linear systems and matrix function approximation in the vein of [34].

Finally, our approach potentially gives rise to new techniques for computing Riemann theta functions. The known asymptotic formulae for the recurrence coefficients of the orthogonal polynomials that we consider consist of the ratio of theta functions composed with the Abel map [12, 40]. In our method, the asymptotics of the recurrence coefficients are essentially governed by the auxiliary function $\mathfrak{h}_{n}$. Since we can efficiently evaluate $\mathfrak{h}_{n}$ in the complex plane, this could potentially enable the efficient evaluation of theta functions. We will also explore this further in future work.