Chebyshev polynomial approximation, multivariate analytic functions, exponential convergence, pluripotential theory, Bernstein–Walsh theorem

41A10, 41A63, 41A25

1 Introduction↩︎

Chebyshev polynomial approximations play a crucial role in many areas of scientific computing, including spectral methods for solving PDEs [1]–[3], the Chebfun software for numerical computing [4], option pricing and finance modeling [5], [6], deep neural networks [7], etc. Consider a function \(f(\mathbf{x})\) defined on the \(d\)-dimensional hypercube \([-1,1]^d\) for some \(d\in\mathbb{N}\), and let \(T_k(x) = \cos(k\arccos(x))\) denote the Chebyshev polynomial of the first kind of degree \(k\). If \(f\) satisfies the Dini–Lipschitz condition [8], then it has a uniformly and absolutely Chebyshev expansion of the form \[\label{onepone} f(\mathbf{x}) = \sum_{\mathbf{k}\in\mathbb{N}_0^d}a_\mathbf{k}T_\mathbf{k}(\mathbf{x}) = \sum_{k_1=0}^{\infty}\cdots\sum_{k_d=0}^{\infty}a_{k_1,\dots,k_d} T_{k_1}(x_1)\cdots T_{k_d}(x_d),\tag{1}\] where \(\mathbb{N}_0^d\) denotes the set of all \(d\)-tuples of nonnegative integers. Truncating this expansion yields the Chebyshev projection of \(f\): \[S_{\Lambda}(f)(\mathbf{x}) = \sum_{\mathbf{k}\in \Lambda}a_{\mathbf{k}}T_{\mathbf{k}}(\mathbf{x}), \notag\] where \(\Lambda \subset \mathbb{N}_0^d\) is a finite multi-index set chosen suitably. In the univariate case with \(\Lambda = \{0,\dots,n\}\), the exponential convergence of the Chebyshev projection for analytic functions was first established by Bernstein in 1912: if \(f\) is analytic in a neighbourhood of \(\overline{E(\rho)}\), where \(E(\rho)\) denotes the open region bounded by the Bernstein ellipse \[\label{oneptwo} \partial E(\rho) := \left\{ z\in\mathbb{C}: z=\frac{u+u^{-1}}{2},\;|u|=\rho \right\}, \quad \rho > 1, \notag\tag{2}\] then \(\max_{x \in [-1,1]} \left|f(x)-S_{\Lambda}(f)(x)\right| = \mathcal{O}(\rho^{-n})\) [9]. The parameter \(\rho\) can be explicitly determined from the locations of the singularities of \(f\). In the multivariate case, however, only a few discussions have been devoted to analyzing the exponential convergence of Chebyshev projection. Key contributions can be summarized as follows:

• Bochner and Martin in [10] established a multivariate analogue of Bernstein’s result, showing that if \(f\) is analytic in a neighbourhood of \(\overline{E(\rho_1)}\times \cdots \times \overline{E(\rho_d)}\), then its Chebyshev coefficients satisfy \(a_{\mathbf{k}} = \mathcal{O}(\rho_1^{-k_1} \cdots \rho_d^{-k_d})\) for \(\mathbf{k} = (k_1, \ldots, k_d) \in \mathbb{N}_0^d\). Although the exponential convergence behavior of Chebyshev projection can be easily seen, the quantification of the exponential convergence rate remains open.

• More recently, based on the observation that multivariate polynomials of fixed degree have anisotropic resolution power in the hypercube, Trefethen in [11], [12] introduced the Euclidean degree for multivariate polynomials, and quantified explicitly the exponential convergence rates of Chebyshev projections with total, Euclidean and maximal degrees. Although the exponential convergence rates for Euclidean and maximal degrees are sharp for Runge functions, e.g., \(f(\mathbf{x})=(x_1^2+ \cdots + x_d^2+h^2)^{-1}\) with \(h > 0\), the sharpness for other types of analytic functions is still unclear (see [13], [14]). Moreover, we will show that even for Runge functions, the exponential convergence rate for total degree is suboptimal (see Remark 6).

The aim of this paper is to analyze the exponential convergence rate of Chebyshev approximation for multivariate analytic functions from the perspective of pluripotential theory. In recent years, pluripotential theory associated to convex bodies has seen considerable development, particularly in connection with multivariate polynomial approximation (see, e.g., [13], [15], [16]). Based on the analyticity assumption determined by the \(P\)-extremal function of the hypercube, we derive a new upper bound estimate for Chebyshev coefficients and show that the exponential convergence of the Chebyshev projection using a downward closed convex polynomial space, which includes commonly used index sets induced by the 1-, 2- and \(\infty\)-norms as special cases, can be quantified by the maximal convergence number. Our result implies that the Chebyshev projection using a downward closed convex polynomial space achieves the best possible exponential rate, since the maximal convergence number actually characterizes the exponential rate of best polynomial approximation for multivariate analytic functions. We further extend our analysis to several related topics, including tensorized Chebyshev interpolation, tensor product Gauss–Legendre quadrature, Padua interpolation and cubature, and Chebyshev-Galerkin method, and establish their exponential convergence rates for multivariate analytic functions as well. Compared with existing convergence results for tensorized Chebyshev interpolation, our new result is sharper (see Remark 8). For Padua interpolation and cubature, to the best of our knowledge, their exponential convergence rates are established for the first time.

The paper is organized as follows. In section 2, we review some basic results on polynomial approximation for multivariate analytic functions. In section 3, we analyze the exponential convergence rate of Chebyshev projection, with supporting numerical experiments provided in section 4. Finally, section 5 applies our analysis results to the related topics, including tensorized Chebyshev interpolation, tensor product Gauss–Legendre quadrature, Padua interpolation and cubature, and Chebyshev-Galerkin method.

2 Bernstein–Walsh theory associated to convex bodies↩︎

In this section, we briefly introduce the Bernstein–Walsh theory associated to convex bodies. We refer to [13], [16], [17] for more details. Throughout the paper, we use the standard multi-index notation, such as \(\mathbf{z}^\mathbf{k}=\prod_{i=1}^d z_i^{k_i}\) for \(\mathbf{z}=(z_1,\ldots,z_d)\) and \(\mathbf{k}=(k_1,\ldots,k_d)\), and \(\boldsymbol{\alpha} \le \boldsymbol{\beta}\) for \(\boldsymbol{\alpha}=(\alpha_1,\ldots,\alpha_d)\) and \(\boldsymbol{\beta}=(\beta_1,\ldots,\beta_d)\) if \(\alpha_i \le \beta_i\), \(1\le i\le d\). We also use \(\Vert \cdot \Vert_{\ell^q}\), \(1 \le q \le \infty\), to denote the standard \(q\)-norm of vectors in \(\mathbb{R}^d\), and \(\Vert \cdot \Vert_{\infty}\) to denote the maximum norm over \([-1,1]^d\).

2.1 The \(\boldsymbol{P}\)-extremal functions and Bernstein–Walsh theorem↩︎

We begin with a review of the classical extremal functions. Let \(PSH(\mathbb{C}^d)\) denote the set of plurisubharmonic functions on \(\mathbb{C}^d\), and let \(L(\mathbb{C}^d)\) denote the usual Lelong class defined by \[L(\mathbb{C}^d) := \left\{ u\in PSH(\mathbb{C}^d) : u(\mathbf{z}) - \max_{1\le i\le d}\{\log|z_i|\} = \mathcal{O}(1) \textrm{ as } |\mathbf{z}| \rightarrow \infty \right\}, \notag\] where \(|\mathbf{z}|\) is the modulus of \(\mathbf{z} = (z_1, \dots, z_d) \in \mathbb{C}^d\). For a compact set \(K\subset\mathbb{C}^d\), the extremal function (or the pluricomplex Green function) of \(K\) is given by \[V_{K}(\mathbf{z}) := \sup \{ u(\mathbf{z}):u\in L(\mathbb{C}^d), \, u \le 0 \textrm{ on } K\}. \notag\] Denote by \(V_K^{*}(\mathbf{z}):=\limsup_{\boldsymbol{\zeta}\rightarrow{\mathbf{z}}}V_K(\boldsymbol{\zeta})\) the upper semicontinuous regularization of \(V_K\). A compact set \(K\) is said to be \(L\)-regular if \(V_{K}\) is continuous on \(K\). Since \(V_K\) is automatically lower semicontinuous, this condition implies that \(V_K^* = V_K\). It is known that the extremal function can be characterized using the subclass of \(L(\mathbb{C}^d)\) arising from polynomials. Specifically, \[V_{K}(\mathbf{z}) = \max \left\{ 0, \, \sup \left\{ \frac{\log|p(\mathbf{z})|}{\textrm{deg}\,p}:\|p\|_{K}:=\sup_{\mathbf{z}\in K}|p(\mathbf{z})| \le 1, \, \textrm{deg}\,p \ge 1 \right\} \right\}, \notag\] where \(p(\mathbf{z})\) is a complex polynomial with total degree \(\textrm{deg}\,p\).

We now introduce the \(P\)-extremal function as a generalization of the classical extremal function. Let \(\mathbb{R}^+=[0,\infty)\) and let \(\Sigma=\{ \mathbf{x} = (x_1,\dots,x_d) \in(\mathbb{R}^+)^d: x_1 + \cdots +x_d \le 1 \}\) be the unit simplex. Consider a convex body (i.e., a compact convex set with nonempty interior) \(P\subset(\mathbb{R}^+)^d\) with the property that \(\Sigma\subset\alpha P\) for some \(\alpha>0\). We define the logarithmic indicator function of \(P\) on \(\mathbb{C}^d\) by \[H_{P}(\mathbf{z}) := \sup_{J\in P\,} \log |\mathbf{z}^{J}| = \sup_{(j_1,\dots,j_d)\in P} \left(j_1\log|z_1|+\cdots+j_d\log|z_d|\right), \notag\] and introduce the generalization of Lelong class as \[\label{Lp} L_P(\mathbb{C}^d) := \left\{ u\in PSH(\mathbb{C}^d) : u(\mathbf{z})-H_P(\mathbf{z}) = \mathcal{O}(1) \textrm{ as } |\mathbf{z}|\!\rightarrow\!\infty \right\}. \notag\tag{3}\] The \(P\)-extremal function of a compact set \(K \subset \mathbb{C}^d\) is then defined as \[V_{P,K}(\mathbf{z}) := \sup \{ u(\mathbf{z}):u\in L_P(\mathbb{C}^d), \, u \le 0 \textrm{ on } K\}. \notag\] Let \(V_{P,K}^*(\mathbf{z}) = \limsup_{\boldsymbol{\zeta}\rightarrow{\mathbf{z}}} V_{P,K}(\boldsymbol{\zeta})\). As in the classical case, we say that a compact set \(K\) is \(PL\)-regular if \(V_{P,K}\) is continuous on \(K\), i.e., \(V_{P,K}^{*}=V_{P,K}\). Note that if \(K\) is \(L\)-regular, then it is \(PL\)-regular for any convex body \(P\) as above. Furthermore, when \(P=\Sigma\), the generalized Lelong class reduces to the usual Lelong class, i.e., \(L_{\Sigma}(\mathbb{C}^d)=L(\mathbb{C}^d)\), and the \(P\)-extremal function recovers the classical extremal function, i.e., \(V_{\Sigma,K} = V_{K}\).

Define the finite-dimensional multivariate polynomial space \[Poly (nP) := \left\{ p(\mathbf{z})=\sum _{\mathbf{k}\in nP \cap \mathbb{N}_0^d} c_{\mathbf{k}} \mathbf{z}^{ \mathbf{k}},\, c_{\mathbf{k}}\in \mathbb{C} \right\} , \quad n \in \mathbb{N}, \notag\] and \[D_n(f,P,K) := \inf \left\{ \Vert f-p \Vert_{K} :p \in \textrm{Poly}(nP) \right\}. \notag\] Note that \(Poly (nP)\) is the space of polynomials of total degree at most \(n\) when \(P=\Sigma\). The Bernstein–Walsh theorem associated to convex body \(P\) is stated as follows [13].

Bernstein–Walsh theorem gives a qualitative result on the best polynomial approximation for analytic functions. Specifically, it shows that \(D_n(f,P,K)\) will decay at the exponential rate \(\mathcal{O}(R_{\scriptsize max}^{-n})\), where \[R_{\scriptsize max} := \sup \{ R:\Omega(P,K,R)\cap S(f) = \emptyset\} = \inf \{ \exp(V_{P,K}(\mathbf{z})):\mathbf{z}\in S(f) \} \notag\] and \(S(f)\) denotes the set of singularities of \(f\). The number \(R_{\scriptsize max}\) is called the maximal convergence number in [18], and it quantifies the exponential convergence rate of best polynomial approximation. We refer to [19] for the history and proofs of the classical Bernstein–Walsh theorem, and to [13] for the proof of the Bernstein–Walsh theorem associated to convex bodies.

2.2 Downward closed convex bodies, nonnegative polars, and the \(\boldsymbol{P}\)-extremal function of the hypercube↩︎

We begin with the definition of downward closed convex bodies and their nonnegative polars.

Below we give two examples of downward closed convex bodies and their nonnegative polars, followed by an example of a convex body that is not downward closed.

• The nonnegative portion of the \(\ell^{q}\) unit ball \(\mathbb{B}_{q}^{+}=\{ \mathbf{x} \in (\mathbb{R}^+)^d: \Vert \mathbf{x} \Vert_{\ell^q} \le 1 \}\) for \(1\le q\le \infty\) is a downward closed convex body. Moreover, the nonnegative polar of \(\mathbb{B}_{q}^{+}\) is \(\mathbb{B}_{p}^{+}\), where \(p\) is the conjugate exponent of \(q\), i.e., \(1/q + 1/p = 1\).

• Let \(L = [0,r_1]\times \cdots \times [0,r_d]\) and \(T = \big\{ \! \sum_{i=1}^{d} \lambda_i \mathbf{e}_i/r_i: \sum_{i=1}^{d} \lambda_i \le 1, \, \lambda_i \ge 0 \big\}\), where \(r_i>0\) for \(i=1,\ldots,d\), and \(\{\mathbf{e}_i\}_{i=1}^d\) is the standard basis of \(\mathbb{R}^d\). Both \(L\) and \(T\) are downward closed convex bodies, and \(L\) is the nonnegative polar of \(T\) (and vice versa).

• Let \(G = \operatorname{conv} \left\{ \mathbf{0}, \mathbf{e}_1, \ldots, \mathbf{e}_d, 2\mathbf{e}_1+2\mathbf{e}_2 \right\}\), where \(d \ge 2\) and \(\operatorname{conv}\) denotes the convex hull. Then \(G\) is a convex body satisfying \(\Sigma\subset G\), but it is not downward closed. Indeed, \(2\mathbf{e}_1 \le 2\mathbf{e}_1+2\mathbf{e}_2 \in G\), whereas \(2\mathbf{e}_1\notin G\).

For a downward closed convex body \(P \subset (\mathbb{R}^+)^d\), we define its centrally symmetric extension to \(\mathbb{R}^d\) by \[C(P) := \left\{ \mathbf{x} \in \mathbb{R}^d: \mathbf{x}^+ := (|x_1|, \dots, |x_d|) \in P \right\}. \notag\] The set \(C(P)\) is convex. Indeed, if \(\boldsymbol{\mu},\boldsymbol{\nu}\in C(P)\), then by the triangle inequality and the convexity of \(P\), we have \[\left(\theta\boldsymbol{\mu}+(1-\theta)\boldsymbol{\nu}\right)^+ \le \theta\boldsymbol{\mu}^+ +(1-\theta)\boldsymbol{\nu}^+ \in P, \quad 0 < \theta < 1,\] which, by the downward closed property of \(P\), implies \(\theta\boldsymbol{\mu}+(1-\theta)\boldsymbol{\nu}\in C(P)\). The gauge function (or the Minkowski functional) of \(C(P)\) is defined as \[\label{norm} \Vert \mathbf{x} \Vert_{P} := \inf \left\{ \lambda\geq0: \mathbf{x} \in \lambda C(P) \right\}, \quad \mathbf{x}\in\mathbb{R}^{d}.\tag{4}\] Since \(C(P)\) is a convex body that is centrally symmetric with respect to the origin, the associated gauge function indeed defines a norm on \(\mathbb{R}^d\). Moreover, it is easily verified that \(P^*\) is also a downward closed convex body. Consequently, we may define the gauge function \(\Vert \cdot \Vert_{P^*}\) in a manner analogous to 4 . The following lemma will be useful.

The Bernstein–Walsh theorem holds for a general compact set \(K \subset \mathbb{C}^d\), yet the closed form of the \(P\)-extremal function \(V_{P,K}\) remains elusive in most cases. However, when \(P\) is a downward closed convex body and \(K\) is a Cartesian product of the form \(K = K_1 \times \cdots \times K_d\), where each \(K_i \subset \mathbb{C}\) is compact and non-polar (i.e., \(K_i\) has a positive logarithmic capacity), it has been shown in [13] that the \(P\)-extremal function of \(K\) can be expressed as \[V_{P,K}^{*}(z_1, \dots, z_d) = \sup_{(y_1,\ldots,y_d)\in P} \left( y_1 V_{K_1}^{*}(z_1) + \cdots + y_d V_{K_d}^{*}(z_d) \right). \notag\] In the case of the hypercube, i.e., \(K_i=[-1,1]\) for \(i=1,\ldots,d\), the extremal function for each \(K_i\) is given by \(V_{K_i}(z)=\log\rho(z)\), in which \(\rho(z) = |z + \sqrt{z^2 - 1}|\) and the branch is chosen so that \(\rho(z) \ge 1\). Consequently, the \(P\)-extremal function of the hypercube \(K = [-1,1]^d\) has the closed form \[\label{p[-1,1]} V_{P,K}(\mathbf{z}) = \sup_{\mathbf{y} \in P} \big(\log {\rho}(\mathbf{z}) \cdot \mathbf{y} \big) = \left\Vert \log {\rho}(\mathbf{z}) \right\Vert_{P^*}, \notag\tag{5}\] where \(\log {\rho}(\mathbf{z}) = \left( \log \rho(z_1), \dots, \log \rho(z_d) \right)\), and the second equality follows from the item (i) of Lemma 1.

3 Exponential convergence analysis of multivariate Chebyshev projection↩︎

In this section, we present a convergence analysis for multivariate Chebyshev projection. We restrict our attention to index sets of the form \(nP \cap \mathbb{N}_0^d\), where \(n \in \mathbb{N}\) and \(P\) is a downward closed convex body. For \(R>1\), we introduce the analyticity domain \[\label{def:AnalRegion} \Omega(P,R):= \{ \mathbf{z} \in \mathbb{C}^d :V_{P,[-1,1]^d}(\mathbf{z}) = \Vert \log \rho (\mathbf{z}) \Vert_{P^*} < \log R \}. \notag\tag{6}\] Note that when \(d = 1\) and \(P = [0,1]\), \(\Omega(P,R)\) reduces to \(E(R)\). In comparison with the two existing analyticity assumptions, namely the Bernstein polyellipse introduced in [10] and the Newton ellipse introduced in [12], our analyticity assumption can lead to the optimal exponential convergence rate.

3.1 A new upper bound of Chebyshev coefficients↩︎

Recall that the Chebyshev coefficient \(a_{\mathbf{k}}\) in (1 ) is defined as \[a_{\mathbf{k}}=\frac{2^{d-\aleph(\mathbf{k})}}{\pi^d} \int_{[-1,1]^d}f(\mathbf{x})T_{\mathbf{k}}(\mathbf{x})\omega(\mathbf{x})\mathrm{d}\mathbf{x}, \notag\] where \(\omega(\mathbf{x})=\prod_{i=1}^d (1-x_i^2)^{-1/2}\) and \(\aleph(\mathbf{k})=\#\{ i:k_i\!=\!0 \}\). The following estimate of Chebyshev coefficients will be used to establish our results, which can be found in [10] and [22]. Here we adopt the latter statement.

The main result of this subsection is stated as follows.

3.2 Exponential convergence rate of multivariate Chebyshev projection↩︎

The following lemma will be used in the convergence analysis for Chebyshev projection.

We now state the main result of this work. Let \(S_{nP}(f)\) denote the multivariate Chebyshev projection of \(f\) with the index set \(nP\cap \mathbb{N}_0^d\), i.e., \[S_{nP}(f)(\mathbf{x})=\sum_{\Vert \mathbf{k} \Vert_{P} \le n}a_{\mathbf{k}}T_{\mathbf{k}}(\mathbf{x}). \notag\] In the following, we establish a quantitative result on the exponential convergence rate of \(S_{nP}(f)\) in the maximum norm over the hypercube \([-1,1]^d\). The notation \(\mathbf{1}\) refers to \((1,\dots,1)\) in the rest of the paper.

Combining Theorem 4 with the Bernstein–Walsh theorem, we conclude that when the index set is induced by a downward closed convex body \(P\), the maximum errors of both multivariate Chebyshev projection and best polynomial approximation converge at the same exponential rate \(\mathcal{O}(R_P^{-n})\). Here, the maximal convergence number \(R_P\) is given by \[\label{RP} R_P := \inf _{\mathbf{z} \in S(f)} \exp(\Vert \log \rho (\mathbf{z}) \Vert_{P^*}),\tag{11}\] where \(S(f)\) denotes the set of singularities of \(f\). Note that the comparison of the convergence rates of Chebyshev, Legendre and their best counterparts in the univariate case (i.e., \(d=1\)) has been comprehensively discussed in [23].

In the special cases where the index sets are induced by 1-, 2-, and \(\infty\)-norm, i.e., \(P=\Sigma\), \(P=\mathbb{B}_{2}^{+}\) and \(P=\Sigma^*\), which correspond to total, Euclidean and maximal degrees respectively, the maximal convergence numbers are respectively given by \[\begin{align} & R_{\Sigma} = \inf _{ \mathbf{z} \in S(f)} \exp\left(\Vert \log \rho (\mathbf{z}) \Vert_{\ell^{\infty}}\right) = \inf _{ \mathbf{z} \in S(f)} \max_{1\leq i\leq d} \{ \rho(z_i) \}, \notag \\ & R_{\mathbb{B}_{2}^{+}} = \inf _{ \mathbf{z} \in S(f)} \exp\left(\Vert \log \rho (\mathbf{z}) \Vert_{\ell^{2}}\right) = \inf _{ \mathbf{z} \in S(f)} \exp\left[ \bigg(\sum_{i=1}^{d} \left( \log\rho(z_i) \right)^2\bigg)^{1/2} \right] , \notag \\ & R_{\Sigma^*} = \inf _{ \mathbf{z} \in S(f)} \exp\left(\Vert \log \rho (\mathbf{z}) \Vert_{\ell^1}\right) = \inf _{ \mathbf{z} \in S(f)} \prod _{i=1}^{d}\rho(z_i). \notag \end{align}\]

4 Numerical experiments↩︎

From Theorem 4 we can deduce that the maximum error of Chebyshev projection of total degree \(n\), i.e., \(S_{n\Sigma}(f)\), decays like \(\mathcal{O}(R_\Sigma^{-n})\), and the maximum error of Chebyshev projection of maximal degree \(n\), i.e., \(S_{n\Sigma^*}(f)\), decays like \(\mathcal{O}(R_{\Sigma^{*}}^{-n})\). In this section, we present several numerical experiments to verify our theoretical results. The maximum errors are estimated as the largest pointwise errors on a uniform tensor-product grid with \(1001\) points in each coordinate direction over \([-1,1]^d\).

5 Applications↩︎

In this section, we apply our analysis results to several closely related topics, including tensorized Chebyshev interpolation, tensor product Gauss–Legendre quadrature, Padua interpolation and cubature, and Chebyshev-Galerkin method. We will show that the exponential convergence rates of these numerical algorithms are determined by the corresponding maximal convergence numbers.

5.1 Tensorized Chebyshev interpolation↩︎

In this subsection, we consider the error analysis of tensorized Chebyshev interpolation. For \(\mathcal{N}=(N_1,\dots,N_d) \in \mathbb{N}_{0}^d\), let \(I_{\mathcal{N}}\) be the polynomial interpolation operator corresponding to the tensor product Chebyshev grid \(\mathcal{S}_1 \times \cdots \times \mathcal{S}_d\), where \(\mathcal{S}_k = \{ \cos (j\pi /N_k) \}_{j=0}^{N_k}\). It is known that \[I_{\mathcal{N}} = I_{x_1}^{N_1} \circ \cdots \circ I_{x_d}^{N_d}, \notag\] where \(I_{x_k}^{N_k}\) denotes the interpolation operator for the variable \(x_k\) at the nodes \(\mathcal{S}_k\). Our error analysis relies on the following lemma, which provides a way to estimate the interpolation error from Chebyshev coefficients.

The following theorem establishes the exponential convergence rate of \(I_{n\mathcal{N}}(f)\) with \(n \in \mathbb{N}\) for multivariate analytic functions.

Theorem 7 shows that the convergence rate of Chebyshev interpolation \(I_{n\mathcal{N}}(f)\) with the tensor product Chebyshev grid \(\mathcal{S}_1 \times \cdots \times \mathcal{S}_d\) is the same as that of Chebyshev projection \(S_{nP}(f)\) with \(P=[0,N_1]\times \cdots \times [0,N_d]\). Both converge at the exponential rate \(\mathcal{O}(R_P^{-n})\) with maximal convergence number \[\label{RN} R_P = \inf _{ \mathbf{z} \in S(f)} \exp\left(\Vert \log \rho (\mathbf{z}) \Vert_{P^*} \right) = \inf _{ \mathbf{z} \in S(f)} \prod_{i=1}^{d} \rho(z_i)^{N_i}.\tag{13}\]

The left subplot of Figure 3 presents the maximum errors of \(I_{n}(f):=I_{n\mathbf{1}}(f)\) for \[\label{testfunc95tencheb} f(x_1,x_2) = \frac{1}{x_1^2x_2^3+1.15}, \quad \frac{1}{x_1^2x_2^4-1.1}, \quad \frac{1}{x_1+x_2+2.1}, \quad \frac{1}{x_1^4+x_2^4-2.05},\tag{14}\] which belong to \(\mathcal{A}_1\), \(\mathcal{A}_2\), \(\mathcal{B}_1\) and \(\mathcal{B}_2\), respectively. Based on \(R_{\Sigma^*}\) in (?? ) and (?? ), the corresponding \(R_{\Sigma^*}\) for these functions are \(\rho(1.15^{1/3})\), \(\rho(1.1^{1/4})\), \(\rho(1.1)\) and \(\rho(1.05^{1/4})\), respectively. We see that \(\Vert f - I_{n}(f) \Vert_{\infty}\) decays exponentially like \(\mathcal{O}(R_{\Sigma^{*}}^{-n})\) for each function, as predicted by Theorem 7.

5.2 Tensor product Gauss–Legendre quadrature↩︎

For \(\mathcal{N}=(N_1,\dots,N_d) \in \mathbb{N}^d\), let \(Q_{\mathcal{N}}(f)\) denote the tensor product Gauss–Legendre quadrature with \(N_i+1\) nodes for the \(i\)-th variable of \(f\), i.e., \[Q_\mathcal{N}(f) = \sum_{k_1=0}^{N_1}\cdots \sum_{k_d=0}^{N_d} \left( \prod_{j=1}^{d} w^{N_j}_{k_j} \right) f\big(x^{N_1}_{k_1},\dots,x^{N_d}_{k_d}\big), \notag\] where \(\{x_k^N\}_{k=0}^N\) are the zeros of the Legendre polynomial of degree \(N+1\), and \(\{w_k^N\}_{k=0}^N\) are the corresponding Gauss–Legendre weights on \([-1,1]\). Let \(Q(f)\) denote the integral of \(f\) over \([-1,1]^d\). The following theorem establishes the exponential convergence rate of \((Q-Q_{n\mathcal{N}})(f)\) with \(n \in \mathbb{N}\) for multivariate analytic functions.

Theorem 10 shows that the tensor product Gauss–Legendre quadrature \(Q_{n\mathcal{N}}(f)\) converges at the exponential rate \(\mathcal{O}(R_P^{-2n})\) for multivariate analytic functions, where the maximal convergence number \(R_P\) is defined in (13 ). The right subplot of Figure 3 illustrates the errors of \(Q_{n}(f):=Q_{n\mathbf{1}}(f)\) for the test functions in (14 ). We see that \(\left| (Q-Q_n)(f) \right|\) decays exponentially like \(\mathcal{O}(R_{\Sigma^{*}}^{-2n})\) for each function, as predicted by Theorem 10.

5.3 Padua interpolation and cubature↩︎

Let \(\mathbb{P}^2_n\) denote the space of bivariate polynomials of total degree at most \(n\), i.e., \(\mathbb{P}^2_n = Poly(n\Sigma)\) with \(\Sigma=\{ (x_1,x_2) \in(\mathbb{R}^+)^2: x_1 +x_2 \le 1 \}\). For \(n \in \mathbb{N}\), the set of the Padua points with cardinality \(\textrm{dim}(\mathbb{P}^2_n) = (n+1)(n+2)/2\) on \([-1,1]^2\) is defined as \[\textrm{Pad}_n := \left\{ \gamma_n\left(\frac{k \pi}{n(n+1)}\right), k = 0,\dots, n(n+1) \right\}, \notag\] where \(\gamma_n(t) = (\cos(nt), \cos((n+1)t))\) traces the algebraic curve given by \(T_{n+1}(x) = T_{n}(y)\) [27]. The Padua points are the first known example of optimal points for total degree polynomial interpolation in two variables, i.e., they are unisolvent and have a Lebesgue constant increasing like \(\mathcal{O}(\log^2n)\) [28].

Let \(I_{\textrm{Pad}_n}\) denote the interpolation operator at the Padua points \(\textrm{Pad}_n\), and let \(Q_{\textrm{Pad}_n}(f)\) denote the associated cubature formula, which is called the nontensorial Clenshaw–Curtis cubature, i.e., \[Q_{\textrm{Pad}_n}(f) = \int_{[-1,1]^2} I_{\textrm{Pad}_n}(f)(x_1,x_2) \textrm{d}x_1 \textrm{d}x_2. \notag\] It is known that \(I_{\textrm{Pad}_n}(f)\) and \(Q_{\textrm{Pad}_n}(f)\) are exact for \(f \in \mathbb{P}^2_n\). We refer to [27], [29] for the derivation and properties of Padua interpolation, and to [28] for those of nontensorial Clenshaw–Curtis cubature. Let \(Q(f)\) denote the integral of \(f\) over \([-1,1]^2\). The following theorem establishes the exponential convergence rate of \(I_{\textrm{Pad}_n}(f)\) and \(Q_{\textrm{Pad}_n}(f)\) for bivariate analytic functions.

Theorem 11 shows that the Padua interpolation and the nontensorial Clenshaw–Curtis cubature converge at the exponential rate \(\mathcal{O}(R_{\Sigma}^{-n})\) for bivariate analytic functions. Figure 4 illustrates the errors of \(I_{\textrm{Pad}_n}(f)\) and \(Q_{\textrm{Pad}_n}(f)\) for \[f(x_1,x_2) = \frac{1}{x_1^2x_2^5+3.5}, \quad \frac{1}{x_1^3+ x_2^3-2.35}, \quad \frac{1}{x_1^2+x_2^2-2.5}, \quad \frac{1}{x_1^2+x_2^2+2}, \notag\] which belong to \(\mathcal{A}_1\), \(\mathcal{B}_1\), \(\mathcal{B}_2\) and \(\mathcal{B}_3\), respectively. Based on \(R_{\Sigma}\) in (?? ), (?? ) and (?? ), the corresponding \(R_{\Sigma}\) for these functions are \(\rho(3.5^{1/7})\), \(\rho(1.175^{1/3})\), \(\rho(1.25^{1/2})\) and \(\rho(\textrm{i})\), respectively. We see that the numerical results are in agreement with our predictions for Padua interpolation. For Padua cubature, whose errors were computed using the Multiprecision Computing Toolbox, the convergence rates are slightly faster than the predicted rates \(\mathcal{O}(R_{\Sigma}^{-n})\).

5.4 Chebyshev-Galerkin method↩︎

Consider the elliptic equation \[\label{elliptic} \nu u(\mathbf{x}) -\Delta u(\mathbf{x}) = f(\mathbf{x}),\quad \mathbf{x}\in \mathcal{D} :=(-1,1)^d,\tag{15}\] with homogeneous Dirichlet boundary conditions \(u|_{\partial \mathcal{D}} = 0\) and \(\nu \ge 0\). Let \(H_{0,\omega}^1(\mathcal{D}) = \{ u\in H_{\omega}^1(\mathcal{D}):u|_{\partial \mathcal{D}} = 0 \}\) be equipped with the norm \(\Vert u\Vert_{1,\omega}=(\Vert u \Vert_{\omega}^2+\Vert \nabla u \Vert_{\omega}^2 )^{1/2}\), where \(\Vert \cdot\Vert_\omega\) denotes the \(L_{\omega}^2(\mathcal{D})\)-norm with Chebyshev weight function \(\omega(\mathbf{x})=\prod_{i=1}^d \left(1-x_i^2\right)^{-1/2}\). A weak formulation of (15 ) is to find \(u\in H_{0,\omega}^1(\mathcal{D})\) such that \[\label{formulation1} a(u,v):=\nu \int_{\mathcal{D}}uv\omega + \int_{\mathcal{D}} \nabla u \cdot \nabla (v\omega) = \int_{\mathcal{D}}fv\omega := (f,v)_{\omega},\quad \forall v \in H_{0,\omega}^1(\mathcal{D}).\tag{16}\] From [3] and the Lax–Milgram lemma [3], it follows that the bilinear form \(a(\cdot,\cdot)\) is continuous and coercive, and the solution of (16 ) exists and is unique. Let \(X_N^0 = \mathrm{span} \big\{ \! \prod_{i=1}^d \psi_{k_i}(x_i):0 \le k_i \le N-2, \, \psi_k(x) = T_k(x)-T_{k+2}(x) \big\}\). The Chebyshev-Galerkin approximation for (16 ) is to find \(u_N \in X_N^0\) such that \[\label{formulation2} a(u_N,v_N) = (f,v_N)_{\omega}, \quad \forall v_N \in X_N^0.\tag{17}\] Using Theorem 4 and some standard arguments, we derive the following convergence result for the numerical scheme (17 ).

To validate our theoretical results, we employ the Chebyshev-Galerkin method proposed in [31] to solve the equation (15 ) for \(d=2\) and \(\nu = 1\), taking \(\widetilde{u}(x,y)\) as the exact solution. Here, \(\widetilde{u}(x,y)\) denotes \(u(x,y)\) with the boundary data lifted according to [32], i.e., \(\widetilde{u}|_{\partial\mathcal{D}} = 0\). It is easy to verify that the maximal convergence number \(R_{\Sigma^*}\) for \(\widetilde{u}(x,y)\) is equal to \(R_{\Sigma^*}\) for \(u(x,y)\). Figure 5 illustrates the \(H_{\omega}^1(\mathcal{D})\)-errors \(\Vert \widetilde{u}-\widetilde{u}_N \Vert_{1,\omega}\) for \[u(x,y) = \frac{\sin(x+y)}{2+x^3 y^4},\quad \frac{4x^2+3y^3}{2.5+x+y},\quad \frac{e^{\sin(x+1.2y)}}{2.35+x^3+y^3},\quad \frac{e^{0.3x+\cos{y}}}{0.09+x^2+y^2}. \notag\] It follows from (?? ), (?? ) and (?? ) that the values of \(R_{\Sigma^*}\) for each \(u(x,y)\) are \(\rho(2^{1/4})\), \(\rho(1.5)\), \(\rho(1.35^{1/3})\) and \(\rho(0.3\mathrm{i})\), respectively. We see that \(\Vert \widetilde{u}-\widetilde{u}_N \Vert_{1,\omega}\) decays exponentially like \(\mathcal{O}(R_{\Sigma^*}^{-N})\) for each case, consistent with our theoretical results.

6 Proofs of (?? ), (?? ) and (?? )↩︎

The idea for the proofs of (?? ), (?? ) and (?? ) is to first obtain a lower bound for \(R_\Sigma\) or \(R_{\Sigma^{*}}\), and then prove that this lower bound can be attained. The following results about \(\rho(z)\) are required.

Acknowledgement↩︎

This work was partially supported by the National Natural Science Foundation of China under grant number 12371367 and the Hubei Provincial Natural Science Foundation of China under grant number 2023AFA083 and the fundamental research funds for the central universities under grant number 2025BRSXA003.

References↩︎