Based on structured matrix techniques, we propose new numerical algorithms to generate a sequence \(\{S_0(x),S_1(x),\dots , S_{K-1}(x) \}\) of Sobolev orthogonal polynomials (SOPs). In this sequence, the polynomial \(S_k(x)\) of degree \(k\) satisfies orthogonality conditions with respect to a Sobolev inner product \(( .,. )\), namely

\[ ( S_k,S_\ell ) \left \{\begin {array}{c} = 0 \textrm { if } k\neq \ell ,\\ \neq 0 \textrm { if } k = \ell . \end {array} \right .\]

A Sobolev inner product is a linear functional where the functions themselves and their derivatives, up to order \(M\), appear with (possibly different) weight functions \(w_m(x)\), i.e.,

\[ (p,q) = \sum _{m=0}^{M} \int _{\Omega } p^{(m)}(x) q^{(m)}(x) w_m(x) dx. \]

Sobolev orthogonal polynomials were first studied in approximation theory, when solving interpolation problems where values of the function and its derivatives are known [1]. The study of their analytical properties, e.g., the behavior of zeros, is an active area of research [4]. Moreover, since the weak form of differential equations gives rise to a Sobolev inner product, SOPs can be used to develop spectral methods [2]. The generation of sequences of SOPs is central to these applications.

Main problem:
Given a Sobolev inner product \((.,.)\), generate the sequence \(\{S_0(x),S_1(x),\dots , S_{K-1}(x) \}\) such that
• \(S_k(x)\) is a polynomial of exact degree \(k\),
• \(( S_k,S_\ell ) = 0\), for \(k\neq \ell \), and \(( S_k,S_\ell ) \neq 0\), for \(k=\ell \).

Our proposed algorithms can be used for general Sobolev inner products, in this presentation we focus on a particular, interesting family of SOPs. Gegenbauer-Sobolev polynomials are orthogonal with respect to the continuous Sobolev inner product

\[ (p,q )_{\mu } = \int _{-1}^{1} p(x) q(x) (1-x^2)^\mu dx + \lambda \int _{-1}^{1} p'(x) q'(x) (1-x^2)^\mu dx. \]

A finite sequence of SOPs \(\{S_0(x),S_1(x),\dots ,S_{K-1}(x)\}\) is also orthogonal to a discrete Sobolev inner product. For Gegenbauer-Sobolev polynomials this discrete inner product can be obtained by applying the Gauss-Jacobi quadrature rule with weights \(\{\alpha _n\}_{n=1}^K\) and nodes \(\{x_n\}_{n=1}^K\) to \((.,.)_\mu \):

\[ (S_k,S_\ell )_{\mu } \approx \langle S_k,S_\ell \rangle _{\mu } = \sum _{n=1}^K \alpha _{n} S_k(x_n) S_\ell (x_n) + \lambda \sum _{n=1}^K \alpha _{n} S_k'(x_n) S_\ell '(x_n). \]

For \(\langle .,.\rangle _\mu \), a Hessenberg matrix \(H_K\) of size \(K\times K\) represents the recurrence relation of the SOPs,

\[ x \left [\begin {array}{ccccc} S_0(x) & S_1(x) & S_2(x) & \dots &S_{K-1}(x) \end {array}\right ] = \left [\begin {array}{cccccc} S_0(x) & S_1(x) & S_2(x) & \dots &S_{K-1}(x) \end {array}\right ] H_K. \]

We show that the matrix \(H_K\) can be computed by solving the following inverse eigenvalue problem [5], where \(H_K\) appears as the \(K\times K\) leading principal submatrix of the solution \(\tilde {H}\).

Inverse eigenvalue problem:
Given a discrete Sobolev inner product \(\langle .,. \rangle \), construct the \(2K\times 2K\) matrices \(\tilde {H}\), \(\tilde {Q}\) such that
• \(\tilde {H}\) is a Hessenberg matrix and \(\tilde {Q}\) is unitary, \(\tilde {Q}^\ast \tilde {Q} = I \).
• The first entries of the unitary matrix \(\tilde {Q}\) are determined by the weights \(\{\alpha _n\}_{n=1}^K\) in the discrete inner product, i.e., \(\tilde {Q} e_1 = w/\Vert w\Vert _2\), where the vector \(w\) contains the weights \(\alpha _n\).
• The matrix \(\tilde {H}\) satisfies the decomposition \(\tilde {H} = \tilde {Q}^\ast J \tilde {Q}\), for a Jordan-like matrix \(J\) determined by the nodes \(\{x_n\}_{n=1}^K\) and their multiplicity, as they appear in the discrete inner product.

For the Gegenbauer-Sobolev polynomials, the vector of weights \(w\) and Jordan-like matrix \(J\) are

\[ w = \left [\begin {array}{c} 0\\ \sqrt {\alpha _1}\\ \vdots \\ 0\\ \sqrt {\alpha _K}\\ \end {array}\right ] \qquad \textrm {and} \qquad J = \left [\begin {array}{ccccccc} x_1 & \sqrt {\lambda }\\ & x_1\\ & & \ddots \\ & & & x_K& \sqrt {\lambda }\\ & & & & x_K\\ \end {array}\right ]. \]

Thanks to this reformulation as a Hessenberg inverse eigenvalue problem, we can use structured matrix techniques to compute \(H_K\). Our proposed numerical algorithm is based on plane rotations and constructs \(H_K\) by employing unitary similarity transformations to the Jordan-like matrix \(J\). It does not require the storage of the whole matrix \(Q_k\), only its first column and is, therefore, more efficient than the state-of-the-art algorithms [3]. Under mild assumptions on the Sobolev inner product, our numerical algorithm can be applied to any sequence of SOPs.

For certain families of SOPs, specialized algorithms can be developed that exploit their properties. For Gegenbauer-Sobolev polynomials we exploit the fact that the Gegenbauer measure forms a symmetrical coherent pair with itself and use properties of (classical) Gegenbauer polynomials to obtain a factorization of \(H_K\) as the product of three structured matrices. The entries of all three matrices are given analytically by closed form expressions, which could reduce the accumulation of rounding error in \(H_K\). We discuss numerical experiments in which we compare the general purpose, specialized, and state-of-the-art algorithms [3] for Gegenbauer-Sobolev polynomials.

References

• [1]  P. Althammer, Eine Erweiterung des Orthogonalitätsbegriffes bei Polynomen und deren Anwendung auf die beste Approximation, J. Reine Angew. Math., 211 (1962), pp. 192–204.

• [2]  L. Fernández, F. Marcellán, T. E. Pérez, and M. A. Piñar, Sobolev orthogonal polynomials and spectral methods in boundary value problems, Appl. Numer. Math., 200 (2024), pp. 254–272.

• [3]  W. Gautschi and M. Zhang, Computing orthogonal polynomials in Sobolev spaces, Numer. Math., 71 (1995), pp. 159–183.

• [4]  F. Marcellán and Y. Xu, On Sobolev orthogonal polynomials, Expo. Math., 33 (2015), pp. 308–352.

• [5]  N. Van Buggenhout, On generating Sobolev orthogonal polynomials, Numer. Math., 155 (2023), pp. 415–443.