Orthogonal polynomials on planar cubic curves

Orthogonal polynomials in two variables on cubic curves are considered, including the case of elliptic curves. For an integral with respect to an appropriate weight function defined on a cubic curve, an explicit basis of orthogonal polynomials is constructed in terms of two families of orthogonal polynomials in one variable. We show that these orthogonal polynomials can be used to approximate functions with cubic and square root singularities, and demonstrate their usage for solving differential equations with singular solutions.


Introduction
We study orthogonal polynomials of two variables with respect to an inner product defined on a planar cubic curve. This is a continuation of recent work by the last two authors that studies orthogonal polynomials on simple one-dimensional geometries embedded in two-dimensional space, including wedges [15] and quadratic curves [16], as well as higher-dimensional cases constructed via surfaces of revolution [17,21,22].
It is assumed the cubic curve γ is of the standard form y 2 = φ(x), where φ is a cubic polynomial of one variable, which includes the standard elliptic curves as a special case. We consider polynomials that are orthogonal with respect to the inner product f, g γ = Ωγ f (x, y)g(x, y)w(x)dσ(x, y), where Ω γ is the set on which the cubic curve γ is defined and w is an appropriate weight function, and the inner product is well defined on the space R[x, y]/ y 2 − φ(x) . These orthogonal polynomials are algebraic polynomials of two variables, but their structure is determined by the characteristics of the curve. In particular, the dimension of the space of the orthogonal polynomials of degree n is 3 for all n ≥ 3, so that it does not increase with n, as with orthogonal polynomials of two variables on either an algebraic surface or on a domain with non-empty interior [5].
Our main result shows that the orthogonal structure on the cubic curve can be understood, by making use of symmetry, through a mixture of two univariate orthogonal systems. To wit, we are able to construct orthogonal polynomials on the curve explicitly in terms of univariate orthogonal polynomials. Moreover, this structural connection propagates to quadrature rules and polynomial interpolation based on the roots of the orthogonal polynomials. Hence, we have developed a toolbox for the computational and analytical study of functions on cubic curves. We provide two applications to showcase the usage of our results. The first is the approximation of functions with cubic singularities. We compare the results to Hermite-Padé approximation, a common technique for approximating functions with singularities, demonstrating that our approximation converges faster and is more robust to degeneracies. The second application is differential equations, which demonstrates the effectiveness of a spectral collocation method, based on our toolbox, for solving differential equations with singular solutions. The examples include the computation of an elliptic integral that can be expressed in terms of the Legendre's incomplete integral of the first kind.
The paper is organized as follows. The orthogonal structure on the cubic curve is established in the next section, where we clarify the families of cubic curves that we consider, which leads to several distinguished cases, and show how orthogonal polynomials can be constructed explicitly in terms of univariate orthogonal polynomials; the section also contains several families of examples. In the third section we consider quadrature rules on the cubic curve as well as polynomial interpolation based on the nodes of the quadrature rules, both from a theoretical and a computational perspective. The applications are discussed in the fourth section and the final section is on possible future work.
Acknowledgment. The first and second authors were supported by the Leverhulme Trust Research Project Grant RPG-2019-144 "Constructive approximation theory on and inside algebraic curves and surfaces".
We consider the cubic curve γ on the plane defined by the standard form since all irreducible bivariate cubics can be transformed 1 into this form [1,2]. We let γ = {(x, y) : y 2 = φ(x)} be the graph of the curve. Without loss of generality, we assume that a 0 > 0, so that φ(x) > 0 for sufficiently large x > 0. Let which is the set on which the cubic curve is defined. The cubic polynomial φ can have either one real zero or three real zeros, so that Ω γ can be either one interval or the union of two intervals. This leads to three possibilities: (I) Ω γ = (A, ∞): the curve has one component; (II) Ω γ = (A 1 , B 1 ) ∪ (A 2 , ∞) with A 1 < B 1 < A 2 : the curve has two disjoint components; (III) Ω γ = (A, B) ∪ (B, ∞): the curve has two touching components.
In the first case, φ has one real zero A. In the second case, φ has three real zeros A 1 < B 1 < A 2 . In the third case, φ has a real zero at A and a double zero at B. For examples, see Figures 1 and 2 below.
One important family of cubic curves included in our definition is that of elliptic curves. An elliptic curve is a plane curve defined by where a and b are real numbers and the curve has no-cusps, self-intersections, or isolated points. This holds if and only if the discriminant is not equal to zero. The graph has two components if ∆ E > 0 and one component if ∆ E < 0. Two elliptic curves are depicted in Figure 1, the left-hand one has one component, whereas the right-hand one has two components.
We need x 2 + Ax + A 2 + a ≥ 0 for x ≥ A, which holds only if its discriminant A 2 − 4(A 2 + a) < 0 or a > −(3/4)A 2 . Under this condition, we have so that the elliptic curve does indeed have one component. Thus, setting c = a + 3 4 A 2 and b = −A 3 − aA we see that the elliptic curve (2.2) becomes If c > 0, then φ(x) = (x−A)(x+ A 2 ) 2 +c has one real zero, so that the elliptic curve has one component. If c < 0, then the curve has two components. One of the advantages of writing the curve in the form (2.3) is that the roots of φ are explicitly given.
We also consider cubic curves that are not elliptic. For example, we can have cubic curves of the form y 2 = a(x 3 − bx 2 ), which will self intersect when b > 0 and will be treated as two components that touch at one point. One example of such curves is depicted in Figure 2.
As an example of the third case, that of closed cubic curves, we mention tear drop curves defined by For a = 1, this curve is inside the unit circle and is depicted in Figure 2.
Our definition also include the curve y 2 = x 3 , which is the case of a = b = 0 in (2.2), but the curve has a singular point and is not an elliptic curve.

2.2.
Orthogonal polynomials on cubic curves. Let y 2 = φ(x) be a cubic curve and let w be a non-negative weight function defined on Ω γ . We consider orthogonal polynomials of two variables that are orthogonal with respect to an inner product defined, on an appropriate polynomial subspace, by where dσ is the arc length measure on the curve. Depending on the support set of w, the integral domain could be compact in the cases I and II. The bilinear form ·, · γ,w defines an inner product on the space R[x, y]/ y 2 − φ(x) . For n ≥ 3, the monomials of degree exactly n, x k y n−k for 0 ≤ k ≤ n, remain of degree n modulo the ring y 2 − φ(x) only when k = 0, 1, 2. In particular, this shows that B n = {y n , xy n−1 , x 2 y n−2 } is a basis of the space of polynomials of degree exactly n in Let V n := V n (γ, w) be the space of orthogonal polynomials of degree n in two variables with respect to this inner product. Applying the Gram-Schmidt process to the basis B n , for example, inductively on n, we obtain the following proposition: Let p n (w) be a family of univariate orthogonal polynomials with respect to the standard inner product on Ω γ , In particular, p n (φw) denotes orthogonal polynomials with respect to φ(x)w(x) on Ω γ . A basis for V n can be given explicitly in terms of orthogonal polynomials with respect to w and φw. We parametrize the inner product (2.4) as More precisely, the domain Ω γ in the integral should be replaced by supp(w) ⊂ Ω γ , where supp(w) denotes the support set of w. For example, in case I, we could choose w so that it has support set [A, B] for some B ∈ R.
We now define an explicit basis for the space V n (w) of orthogonal polynomials on the cubic curve γ. We denote this basis by Y n,i and denote its squared norm by H n,i . The squared norm of p n (w; x) is denoted by h n (w) = p n , p n w .
Theorem 2.2. Let γ be a cubic curve and let w be a weight function defined on Ω γ .
1. For n = 0 and n = 1, we define Then V 0 = span{Y 0 } and V 1 = span{Y 1,1 , Y 1,2 }. Moreover, 2. For m ∈ N 0 and m ≥ 1, we define Then Y n,i is a polynomial of degree n in R[x, y]/ y 2 − φ(x) for i = 1, 2, 3 and Moreover, the norms of these polynomials are given by Proof. For n ≥ 2, we first show that Y n,i is of degree n in R[x, y]/ y 2 −φ(x) . Throughout this proof, we introduce a function ψ(x, y) so that the equation of the cubic curve becomes Now, for n = 2m, we write which is a polynomial of degree 2m in x, y variables. The same argument also shows that shows that Y 2m,3 is a polynomial of degree 2m mod y 2 − φ(x) . A similar argument works for Y 2m+1,i . We now verify orthogonality. The Y n,i are of two forms: either Y n,i (x, y) = f (x) or Y n,i (x, y) = yg(x). due to symmetry in the expression of ·, · γ,w in (2.6), This establishes orthogonality between all Y n,i of the form f (x) and those of the form yg(x). To prove orthogonality between the remaining Y n,i (those that are both of the form f (x) and those that are both of the form yf (x)), we relate the bivariate inner product on γ to the univariate inner product on Ω γ in both cases: observe that Thus, the orthogonality between the remaining Y n,i follow from the orthogonality of the p n (w) and p n (φw). Hence we have showed that the Y n,i form a basis for V n , n ≥ 0. It follows from the last two equations that if Y n,i = p k (w), then H n,i = Y n,i , Y n,i γ,w = 2 p k (w), p k (w) w = 2h k (w). If Y n,i = yp k (φw), then H n,i = 2 p k (φw), p k (φw) φw = 2h k (φw).
2.3. Fourier orthogonal series. For w defined on R, the Fourier orthogonal series in terms of orthogonal polynomials {p n (w)} is defined by where the identity holds in L 2 (w) as long as polynomials are dense in L 2 (w), which we assume to be the case. Furthermore, let s n (w; f ) denote the n-th orthogonal partial sum of this expansion; that is, Likewise, let γ be a cubic curve and w be a weight function defined on γ, we define the Fourier orthogonal series of f ∈ L 2 (γ, w) by The n-th partial sum of this expansion is denoted by S n (w; f ); that is The next theorem shows that this partial sum can be written in terms of the partial sums of orthogonal series with respect to w and φw. Let · w denote the norm of L 2 (γ, w).
Theorem 2.3. Let γ be a cubic curve and let w be a weight function defined on Ω γ . (2.11) Proof. It follows from the definitions of f e and ·, · γ,w that , we can write the partial sum with n = 2m as A similar proof works for S 2m+1 (w; f ). By (2.7) and the Parseval identity, we see that where the second identity follows since |s 3m (w; f e , x)| 2 does not contain y, whereas |ys 3m (φw; f o , x)| 2 contains a y 2 , which is equal to φ(x). The proof for the norm of S 2m+1 (w; f ) is similar. This completes the proof.
Corollary 2.4. Let γ be a cubic curve and let w be a weight function defined on Ω γ .
Proof. For f ∈ L 2 (γ, w), it follows from (2.6) and |a + b| 2 ≤ 2(|a| 2 + |b| 2 ) that The general theorem of orthogonal polynomials of several variables shows that where A 0,i are 1 × 2 matrices and B n,0 is a real number; A 1,i are 2 × 3 matrices and B 1,i are 2 × 2 matrices; A n,i and B n,i are 3 × 3 matrices for all n ≥ 2. The matrices A n,i and B n,i are determined by orthogonality relations: In particular, by (2.7), (2.8) and (2.9) it is easy to see that these matrices are of the form and These three-term relations in two variables hold when (x, y) are on the cubic curve γ or modulo the polynomial ideal y 2 − φ(x) . It is worth mentioning that we obtain, for example, where (A) i,j stands for (i, j)-element of the matrix A, in which the lefthand side is a polynomial of degree n + 1, where the righthand side is of degree n. This holds without contradiction because (x, y) ∈ γ.

2.5.
Examples of orthogonal polynomials on cubic curves. Recall from Theorem 2.2 that we require two one-variable orthogonal polynomial families to construct an orthogonal basis on the cubic curve. We shall always choose w to be either the Jacobi or Laguerre weight. As we shall argue below the second family of orthogonal polynomials p n (φw) will be non-classical if at least one of the roots of φ is not at an endpoint of supp(w). Otherwise, if all the roots of φ are at the endpoint(s) of supp(w), the second orthogonal polynomial family will also be Jacobi or Laguerre but with different parameters for its weight function. The following cases arise: 1. φ has one simple real root, or φ = (x − a)g 2 (x), where g 2 (x) is a second degree polynomial with a pair of complex conjugate roots. 2. φ has one triple real root, or φ = k(x − a) 3 , where k is a constant. 3. φ has a simple and a double real root, or b < a. The latter case is not of interest since then b represents an isolated point of the curve y 2 = φ. 4. φ has three distinct real roots, In each of these cases the domain Ω γ on which the orthogonal polynomials are defined is either a semi-infinite interval [A, ∞) or the union of a compact interval [a, b] and a semi-infinite interval. On semi-infinite intervals, we additionally consider the cases for which the support of the weight is (i) also semi-infinite, or (ii) compact. Through a linear change of variables we may, without loss of generality, let supp(w) ⊂ Ω γ be the canonical compact interval [−1, 1] or the semi-infinite interval [0, ∞). If Ω γ is the union of a compact and semi-infinite interval, we may consider [−1, 1] and [0, ∞) separately. If supp(w) includes the root(s) of φ, then these will be mapped to −1 and/or 1 if the domain is compact and to 0 if the domain is semi-infinite. This implies that on [−1, 1], if we let the weight w be the Jacobi weight w α,β , then there exists a polynomial g k (x) of degree k, 0 ≤ k ≤ 3, that is strictly positive on supp(w) and whose roots are outside supp(w), such that Hence, φ has i + j roots at the endpoints of [−1, 1] and k roots outside [−1, 1]. The weight function of the second orthogonal polynomial family is φw = w α+i,β+j g k , which is a non-classical weight if k > 0 and classical if k = 0. On [0, ∞), if we let w be the Laguerre weight w α , then Here φ has j roots at the endpoint of [0, ∞) and k roots outside [0, ∞). The weight function of the second family is φw = w α+j g k , which is non-classical if k > 0 and classical if k = 0.
We consider three cubic curves as examples. For the first two, orthogonal polynomials can be given explicitly in terms of classical orthogonal polynomials. The third example discusses orthogonality on elliptic curves.
2.5.1. Orthogonal polynomials on the curve y 2 = x 3 . In this example, the curve and the weight functions are Since all the roots of φ are at 0, the endpoint of supp(w), the orthogonal polynomial basis on the curve can be constructed entirely out of Laguerre polynomials. The polynomial p n (w; is also an Laguerre polynomial with parameter α +3. In this setting the inner product on the curve becomes The orthogonal basis B n of the space V n in Theorem 2.2 becomes The norm of the Laguerre polynomial L (α) n is given by from which the norm of the basis in V n can be derived as in Theorem 2.2.
2.5.2. Jacobi polynomials on tear drop curves. In this example, the curve is the tear drop curve and the weight function is the Jacobi weight, for α, β > −1, Since all the roots of φ are at the endpoints of supp(w α,β ), the orthogonal basis on the tear drop curve can be constructed entirely out of Jacobi polynomials. The polynomial p n (w; x) is the usual Jacobi polynomial P and p n (φw) is also a Jacobi polynomial, p n (φw) = P (α+2,β+1) n . In this setting the inner product on the curve becomes The orthogonal basis B n of the space V n in Theorem 2.2 becomes

2.5.3.
Orthogonal polynomials on elliptic curves. We consider two elliptic curves. The first one is given by which has one component. We set u = x + 2 so that φ = u (u − 3) 2 + 1 has a root at u = 0 and we can choose the classical Laguerre weight w α (u) = u α e −u , α > −1, defined for u ≥ 0. In this setting, the {p n (w)} are given by the Laguerre polynomials L (α) n . Since φ(u) has two roots outside supp(w α ), the orthogonal polynomial family {p n (φw)} is non-classical and orthogonal with respect to The inner product in this setting becomes If we choose a weight function w that is supported on x ∈ [−2, 2], say, then the inner product is defined on a finite segment of γ. We set x = 2u, where u ∈ [−1, 1] and choose the Jacobi weight w = w α,β (u). In this case, the {p n (w)} are given by the Jacobi polynomials P (α,β) n and the {p n (φw)} are non-classical orthogonal polynomials (since φ has two roots outside supp(w)) with respect to the weight φ(u)w α,β (u) = 2 (2u − 1) 2 + 1 w α,β+1 (u).
Our second elliptic curve is given by which has two components. The first one is a closed curved with −2 ≤ x ≤ 0 and the second one is an open curve defined for x ≥ 2. On the first component we set u = x+1 ∈ [−1, 1] and choose w = w α,β (u) in which case p n (w; u) = P (α,β) n (u) and the polynomials {p n (φw; u)} are orthogonal with respect to φw = w α+1,β+1 (u − 3). On the second component we let u = x − 2 and choose w = w α (u) so that p n (w; u) = L (α) n (u) and the polynomials {p n (φw; u)} are orthogonal with respect to φw = w α+1 (u + 4)(u + 2).
Remark: For a particularly convenient and flexible method for numerically computing these non-classical orthogonal polynomials p n (φw), based on the Lanczos algorithm, see [14,13]. This is described in more detail in Section 4.

Quadrature rules and polynomial interpolation
We consider quadrature rules and polynomials interpolation on the cubic curve.
where λ k,N are the Gaussian quadrature weights. Let γ be a cubic curve. We denote by Π n (γ) the space of polynomials of degree at most n restricted to the curve γ. By Theorem 2.2, 3m−2 . From Proposition 2.1 it follows dim Π 0 (γ) = 1 and dim Π n (γ) = 3n, n ≥ 1.
The quadrature (3.2) on the cubic curve is an analogue of the Gaussian quadrature rule on the real line. We now consider polynomial interpolation based on the nodes of this quadrature rule.

Lagrange interpolation.
First we recall the univariate Lagrange interpolation polynomial on the zeros x k,N , 1 ≤ k ≤ N , of p N (w), denoted by L N (w; f ), which is the unique polynomial of degree at most N − 1 that satisfies for any continuous function f . It is well-known that L N (w; f ) is given by .
By the Christoffel-Darboux formula, we can also write k as Furthermore, it is the unique interpolation polynomial in Π n (γ) if n = 2m + 1 and in Proof. If n = 2m + 1, then N = 3m + 1 and we interpolate at 2N = 6m + 2 points. In this case, both L N (w; f e ) and L N (w; f o ) are elements of Π 3m , it follows that L n (w; f ) ∈ Π 2m+1 (γ). Since Y 2m+1,1 (x, y) = p 3m+1 (w; x) vanishes on all interpolation points, there are dim Π 2m+1 (γ) − 1 = 6m + 2 independent functions over the set of nodes in the space Π 2m+1 (γ). Similarly, if n = 2m, then N = 3m and we interpolate at 2N = 6m points. In this case, both L N (w; f e ) and L N (w; f o ) are elements of Π 3m−2 , it follows that L n (w; f ) ∈ Π 2m (γ) ∪ {Y 2m+1,3 } since Y 2m+1,3 = yp 3m−1 (w; x). Since Y 2m,1 (x, y) = p 3m (w; x) vanishes at all nodes, we see that there are dim Π 2m (γ) − 1 + 1 = 6m independent functions over the set of nodes in the space. Now, for 1 ≤ k ≤ N , we obtain from the Lagrange interpolation of L N (w; f ) that N , y k,N ); similarly, we also obtain that and ·, · N denotes the discretised univariate inner product ·, · w based on the N -point Gaussian quadrature rule: where, as before, x k,N and λ k,N are the Gaussian quadrature nodes (the roots of p N (w)) and weights, respectively. According to the following bivariate analogue of the univariate result just mentioned, we require a system of M orthogonal functions with respect to an M -point discrete inner product to construct an interpolant expanded in an orthogonal basis. We first consider the inner product ·, · N,γ coming from discretization of ·, · γ,w via Gauss quadrature,  and with n = 2m and N = N n = 3m, the 2N − 1 functions are the largest sets of functions from among the Y k,i defined in Theorem 2.2 that are orthogonal and have nonzero norms with respect to ·, · N,γ .
Proof. From the symmetry of ·, · N,γ , we have, similar to the property (2.7) of the continuous inner product, This proves the orthogonality between the basis functions Y k,i of the form f (x) and those of the form yg(x). To demonstrate orthogonality between the remaining Y k,i (those that are both of the from f (x) or both of the form yf (x)), we note that since N -point Gaussian quadrature is exact for polynomials of degree ≤ 2N − 1, (3.11) for 0 ≤ k, j ≤ N − 1. Since φ has degree 3, for 0 ≤ k, j ≤ N − 2. Equations (3.11) and (3.12) demonstrate the orthogonality of the sets of functions in (3.10) and (3.9) as well as their nonzero norms. The aforementioned bounds on the indices k and j and the fact that p N (w) vanishes at all the quadrature nodes imply that it is not possible to add a function p k (w), with k ≥ N or yp k (φw), with k ≥ N − 1 to the sets (3.10) or (3.9) such that it is orthogonal and have a nonzero norm with respect to ·, · N,γ . In particular, it is not possible to add yp N −1 to (3.10) or (3.9) We conclude that it is not possible to construct an interpolant via Gauss quadrature in the manner of Proposition 3.3. It is nevertheless possible to construct an interpolant a la Proposition 3.3 if ·, · γ,w is discretised via Gauss-Radau or Gauss-Lobatto quadrature [7].
As shown in the proof of Theorem 3.2, the set of 2N functions {p k (w), yp n (φw)} N −1 k=0 at the 2N nodes (x k,N , ±y k,N ), k = 1, . . . , N are linearly independent. Hence, the interpolant in Theorem 3.2 can be represented in the form (3.13) and the coefficients a k,N and b k,N can be obtained by solving a 2N ×2N Vandermondelike linear system whose columns consist of p k (w) and yp n (φw), evaluated at (x k,N , ±y k,N ) for k = 1, . . . , N . However, since this procedure requires O(N 3 ) operations, we introduce an alternative method with which the coefficients can be obtained with a complexity of either O(N 2 ) (for a general weight w) or O(N log N ) (for the Chebyshev weight).
We now introduce a new inner product and orthogonal basis with respect to which it is possible to construct an interpolant with Gaussian quadrature. Although the inner product ·, · γ,w is the natural parametrization of (2.4), the new inner product gives rise to a new orthogonal basis that is more efficient for computational purposes compared to the basis in Theorem 2.2, as we shall demonstrate.
We define the inner product [·, ·] γ,w as where f e , g e , f o , g o are the even and odd parts, respectively, of functions f (x, y) and g(x, y) on γ, as defined in (2.10). Note that the odd part of functions on γ have removable singularities at points (x, y) where y = φ(x) = 0.
In the inner product space defined by [·, ·] γ,w , the orthogonal polynomial basis on γ is exactly the same as in Theorem 2.2 (where the inner product space is equipped with ·, · γ,w ), except that the Y n,i that are of the form yp k (φw) are replaced by yp k (w). That is, whereas the orthogonal polynomial basis in Theorem 2.2 is constructed from p n (w) and p n (φw), the orthogonal basis with respect to [·, ·] γ,w is constructed solely out of p n (w). Hence, in the inner product space with [·, ·] γ,w , one only requires univariate classical orthogonal polynomials. The results in Theorem 2.3 and Corollary 2.4 also hold for the orthogonal basis with respect to [·, ·] γ,w , mutatis mutandis. A notable difference between the inner products ·, · γ,w and [·, ·] γ,w is that, unlike the Jacobi operators for ·, · γ,w in section 2.4, the latter gives rise to a non-symmetric Jacobi operator for multiplication by y. This is because the inner product [·, ·] γ,w is not self-adjoint with respect to multiplication by y. That is, 3.4. Interpolation via quadrature with respect to [·, ·] γ,w . Discretising [·, ·] γ,w via Gauss quadrature, we obtain We shall need the following result.

15)
and with n = 2m and N = N n = 3m, the 2N functions Proof. It follows from the definition (2.10) that for a function on γ that depends only on x, and for a function on γ of the form yg(x), .
Hence, [f (x), yg(x)] N = 0, which proves orthogonality between functions in the sets (3.16) and (3.15) of the form p k (w) and those of the form yp k (w). We note that since the N -point Gaussian quadrature rule is exact for polynomials of degree ≤ 2N − 1, it follows that a k,N p k (w; x) + y and where N = N n = 3m + 1 if n = 2m + 1 and N = N n = 3m if n = 2m.
Proof. It follows from Proposition 3.5 that for 0 ≤ k ≤ N − 1, By the previously mentioned result for univariate interpolants in section 3.3, this implies that Note that if w is chosen to be the Chebyshev weight, one can compute the coefficients of the interpolant, a k,N and b k,N , in O(N log N ) operations using the fast cosine transform. A fast transform is also available for the uniform Legendre weight w = 1 [10], and potentially for any Jacobi weight [9]. Since fast transforms are not available for the non-classical polynomials p n (φw), O(N 2 ) operations are required to compute the coefficients of the interpolant obtained via (Gauss-Radau or Gauss-Lobatto) quadrature. The inner product [·, ·] γ,w is therefore preferable to ·, · γ,w for computing interpolants via quadrature.

Applications
We can approximate functions of the form which have square root-type singularities at the zero(s) of φ(t), by recasting them as functions f (x, y) on the cubic curve γ = (x, y) : y 2 = φ(x) . If f (x, y) is a smooth function of x and y on γ, then the bivariate interpolant on γ will converge much faster compared to the univariate interpolant of g(t). Similarly, if φ has an inverse φ −1 on an interval, then we can approximate functions of the form which have cubic-type singularities where φ −1 has zeros, by an interpolant f (x, y) on γ by setting y = t and x = φ −1 (t 2 ). First we approximate a function of the form (4.2) and in the next section we consider functions of the form (4.1) that arise as solutions to differential equations.
To compute p k (w) (and p k (φw)) we use a variant of the Stieltjes procedure where orthogonal polynomials are calculated via computing the connection coefficients with Legendre polynomials, that is, it computes the expansion p k (w) = k j=0 c kjPk whereP k are normalized Legendre polynomials. Provided w is a polynomial then the inner products are computable exactly via the Legendre Jacobi operator J, that is, Therefore we can readily determine c kj by orthogonalizing via Gram-Schmidt 2 . This representation makes differentiation straightforward as we know P k (x) = (k+1)/2P (18.9.15)]. In practice, we take the Legendre weight w(x) = 1 so that p k (w) =P k and we need only calculate p k (φ). 2 Equivalently, this can be viewed as Lanczos iteration with a non-standard, banded inner product as explained in [13]. A convenient implementation is available in the Julia package OrthogonalPoly-nomialsQuasi.jl [14] 4.1. Function approximation. To approximate the function where J denotes the Bessel function of the first kind, we can set y 2 = φ(x) = x 3 − 2 , hence φ −1 (y 2 ) = 3 y 2 + 2 . Then For comparison purposes with standard bases, we also approximate (4.3) using algebraic Hermite-Padé (HP) approximation [6]. Given the function values f (x k,N ), p 0 + p 1 f + p 2 f 2 + · · · + p m f m N = minimum. Here, the norm · 2 N := ·, · N is induced by the following discrete inner product We assume some kind of normalization so that the trivial solution p 0 = . . . = p m = 0 is not admissible. The HP approximant of f (x), viz. ψ(x), is the algebraic function defined by (4. 7) p 0 (x) + p 1 (x)ψ(x) + p 2 (x)ψ 2 (x) + · · · + p m (x)ψ m (x) = 0.
In practice, we compute the polynomials p 0 , . . . , p m by expanding them in an orthonormal polynomial basis with respect to the discrete inner product (4.6). Then (4.5) reduces to a least squares problem whose solution we compute with the SVD. Hence the implicit normalization used is that the vector of polynomial coefficients of p 0 , . . . , p m in the orthonormal basis is a unit vector. This computational approach is similar to that used in [8,19] for the case m = 1, which corresponds to rational interpolation or least squares fitting. Throughout we shall consider diagonal HP approximants for which the degrees of the polynomials p 0 , . . . , p m are equal, say degree d. We require that the number of points at which f is sampled is greater than or equals to the number of unknown polynomial coefficients: If N = m(d + 1) + d, then the minimum attained by the solution to the least squares problem (4.5) is zero (4.8) N = m(d + 1) + d, ⇒ p 0 + p 1 f + p 2 f 2 + · · · + p m f m N = minimum = 0. We call this the interpolation case. If N > m(d+1)+d, then the minimum attained by the least squares solution will be nonzero in general. Throughout we shall approximate functions with HP interpolants.
Note that if m = 1 and p 1 (x) = 1 in (4.8), then the HP approximant, ψ(x) = −p 0 (x), is a polynomial interpolant of f on the grid; if m = 1, then the HP approximant, ψ = −p 0 (x)/p 1 (x), is a rational interpolant of f (with poles in the complex x-plane). If m ≥ 2, then for every x, ψ(x) will generally be an m-valued approximant of f (with poles and algebraic branch points in the complex x-plane). We want to pick only one branch of the m-valued function ψ to approximate f . One way to do this is to solve (4.7) with Newton's method using a polynomial or rational approximant as first guess. Figure 3 compares the rate of convergence to g(t) in (4.3) of HP approximants with m = 0, 1, 2, 3 (polynomial, rational, quadratic and cubic HP interpolants) and interpolants on the cubic curve (4.4) which were obtained via quadrature in the Chebyshev basis as described in section 3.4. The figure shows that the interpolant on the cubic curve γ converges super-exponentially (since f is an entire function in x and y), which is significantly faster than HP approximants (which in addition appear to have stability/ill-conditioning issues). Moreover, it is robust as → 0, whereas HP interpolants break down in accuracy.  for y 2 = φ(x), subject to boundary conditions. As discussed in section 2.5, we can let Recall that φ may vanish at the endpoint(s) of Ω γ . If φ vanishes on supp(w), which can only happen at the endpoint(s) of supp(w), the differential equation is singular in x. This is to be expected since linear differential equations with singular solutions must be singular. However, this does not imply that the solution will be singular on γ. For example, the Bessel function (4.3) is the solution to a differential equation of the form (4.9) that is a singular function of x but an entire function of x and y on the curve (4.4).
For the orthogonal basis with respect to the inner product [·, ·] γ,w , which consists of only p n (w), solutions to (4.9) can be computed with the ultraspherical spectral method [12] if the p n (w) are chosen to be the Chebyshev polynomials. This approach represents (4.9) as banded matrices in coefficient space, in contrast to the dense matrices that arise in collocation methods. It is also possible to devise a spectral method with banded operators for the basis consisting of p n (w) and p n (φw). For this basis, unlike the basis consisting solely of p n (w), the operators are not known explicitly but can be constructed efficiently with quadrature.
In the examples that follow, however, we compute solutions to (4.9) with a spectral collocation method since we found that it was just as accurate (for the first example below) or more accurate (for the second example) than the ultraspherical spectral method 3 . In particular, we approximate the solution as u n k,i Y n,i + u n n,2 Y n,2 + u n n,3 Y n,3 , (4. 10) if n is odd (cf. (3.15)) and u(x, y) ≈ u n 0 Y 0 + u n 1,1 Y 1,1 + u n 1,2 Y 1,2 (4.11) if n is even (cf. (3.16)). It follows from (3.15) and (3.16) that an equivalent representation of the right-hand sides of (4.10) and (4.11) is given by the right-hand side of (3.17) (for the basis consisting of only p n (w)) or (3.13) (for the basis constructed from p n (w) and p n (φw)). The coefficients u n k,i (or, equivalently, a k,N and b k,N in (3.17) or (3.13)) are determined by solving a linear system so that the ODE is satisfied on a grid. We let the grid be (x k,N , ±y k,N ), k = 1, . . . , N , where the x k,N are the roots of 3 The well-conditioning of the ultraspherical method is established in [12] for non-singular linear differential equations. Numerical evidence in [3] illustrated superior conditioning and accuracy of the ultraspherical method compared to collocation for a linear equation with regular singularities. An analysis of the ultraspehrical and collocation method for equations with irregular singularities (which can include equations of the form (4.9)) remains an open problem which is beyond the scope of the present paper. p N (w) and y k,N = φ(x k,N ). Recall that if n = 2m + 1, then N = 3m + 1 and if n = 2m, then N = 3m.
To implement the collocation method, we use the fact that and higher order derivatives can be derived recursively in a similar manner. Note that since the collocation points x k,N are in the interior of supp(w), φ(x k,N ) = 0 and hence (4.12) is well defined at the collocation points. subject to u(−1, φ(−1)) = 0 and evaluated at (x, y) = (1, φ(1)). For = 0, the integral can be computed efficiently with Gauss-Jacobi quadrature because the singularity of the integrand at x = −1 is incorporated into the weight and the remainder of the integrand is analytic on [−1, 1]. If the singularity is close to −1, i.e., for 0 < 1, Gauss-Jacobi quadrature is expensive because one needs a high degree polynomial to resolve the nearly singular integrand. By solving the differential equation on the curve y 2 = φ, the computational cost of computing the integral is essentially independent of because the coefficients of the differential equation are polynomials (and hence analytic) on the curve for all ≥ 0. Figure 4 illustrates a solution to (4.13) which is accurate to 15 digits. and satisfies a second order equation of the form (4.9) with yφa 2 (x, y) = c 2 x 2 φφ + c 2 φ 2 + yc 1 φ ∈ Π 6 (γ), ∂ y a 2 (x, y) ∈ Π 5 (γ), yφa 0 (x, y) = yφa 3 2 (x, y) ∈ Π 9 (γ).
(4.15) Figure 5 shows f (x) and the accuracy of the approximate solutions obtained by solving the differential equation, subject to Dirichlet boundary conditions, with a collocation method using two orthogonal bases (one in the inner product space defined by ·, · γ,w and the other in the space defined by [·, ·] γ,w ). The accuracy obtained with both bases are comparable and roughly four digits of accuracy are lost, most likely due to the large condition numbers of the matrices that arise in the collocation methods. The coefficients decay super-exponentially because the solution is an entire function of x and y on the γ. By contrast, a conventional spectral method using an orthogonal basis on [−1, 1] could only achieve algebraic convergence to the solution because of its singularities at x = ±1.

Conclusion
We have constructed orthogonal polynomials on cubic curves and demonstrated that they can be used for function approximation and solving differential equations. There are some clear questions left to explore: (1) Higher order algebraic curves. In the case where there is symmetry upon reflection across the x axis we may be able to reduce to one-dimensional orthogonal polynomials, but a general construction is as-of-yet unclear. (2) Applications to partial differential equations. We mention that the construction of OPs on quadratic curves in [16,17] led to numerical methods for partial differential equations on non-standard disk-slices and trapeziums in [20] and explicit formulae for inverting wave operators on a cone [18]. This suggests that the OPs on cubic curves introduced here may be applicable to solving partial differential equations whose boundaries are defined in terms of a cubic curves. with c 1 = 10 and c 2 = 5, which is the solution to a second order differential equation of the form (4.9) with coefficients given in (4.15). Top-right: the coefficients (defined in (3.17) and (3.13)) obtained by solving the differential equation with a collocation method in the basis consisting of (i) p n (w) and yp n (w) and (ii) p n (w) and yp n (φw) with w = 1. Bottom: the maximum error of the collocation solutions on [−1, 1] with 2N coefficients.