On Sobolev bilinear forms and polynomial solutions of second-order differential equations

Given a linear second-order differential operator L≡ϕD2+ψD\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {L}}\equiv \phi \,D^2+\psi \,D$$\end{document} with non zero polynomial coefficients of degree at most 2, a sequence of real numbers λn\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _n$$\end{document}, n⩾0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n\geqslant 0$$\end{document}, and a Sobolev bilinear form B(p,q)=∑k=0Nuk,p(k)q(k),N⩾0,\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\mathcal {B}}(p,q)\,=\,\sum _{k=0}^N\left\langle {{\mathbf {u}}_k,\,p^{(k)}\,q^{(k)}}\right\rangle , \quad N\geqslant 0, \end{aligned}$$\end{document}where uk\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbf {u}}_k$$\end{document}, 0⩽k⩽N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0\leqslant k \leqslant N$$\end{document}, are linear functionals defined on polynomials, we study the orthogonality of the polynomial solutions of the differential equation L[y]=λny\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {L}}[y]=\lambda _n\,y$$\end{document} with respect to B\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {B}}$$\end{document}. We show that such polynomials are orthogonal with respect to B\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {B}}$$\end{document} if the Pearson equations D(ϕuk)=(ψ+kϕ′)uk\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D(\phi \,{\mathbf {u}}_k)=(\psi +k\,\phi ')\,{\mathbf {u}}_k$$\end{document}, 0⩽k⩽N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0\leqslant k \leqslant N$$\end{document}, are satisfied by the linear functionals in the bilinear form. Moreover, we use our results as a general method to deduce the Sobolev orthogonality for polynomial solutions of differential equations associated with classical orthogonal polynomials with negative integer parameters.


Introduction
For a fixed integer N ≥ 0, let u k , 0 ≤ k ≤ N , be linear functionals defined on polynomials. Then, a symmetric bilinear form B can be defined as (1.1) A bilinear form involving derivatives, such at (1.1), is said to be a Sobolev bilinear form. The case when u N is a regular (or quasi-definite) linear functional and u k , 0 ≤ k ≤ N −1, are Dirac deltas supported on convenient points of the real line has been of great interest. In [9], the following discrete-continuous Sobolev bilinear form was studied: where c 0 , . . . , c N −1 are distinct real numbers, N is a fixed positive integer, u N is a regular linear functional, and A is an N × N quasi-definite symmetric real matrix, that is, its principal submatrices are nonsingular. Notice that (1.2) can be written as where D I and C I are bilinear forms associate with the discrete and continuous parts of B I , respectively. Using the fact that A is quasi-definite, it was shown in [9] that if the monic polynomials {Q n (x)} n≥0 are orthogonal with respect to B I , then, for n ≥ N : n=0 are orthogonal with respect to D I . Similar results were obtained in [10], where the bilinear form (1.1) was considered with u N regular, and taking u k , 0 ≤ k ≤ N − 1, to be derivatives of Dirac deltas supported on some points of the real line.
Special attention has been given to the case when u N is associated with the classical Jacobi or Laguerre orthogonal polynomials, and u k , 0 ≤ k ≤ N − 1, are derivatives of Dirac deltas supported on the endpoints of the interval of orthogonality of u N (see [1][2][3]5,15,18,19] and the references therein). Such bilinear forms arise when studying the orthogonality of Jacobi or Laguerre polynomials with the so-called non standard parameters, that is, negative integer parameters such that the coefficient c n in the corresponding three-term recurrence relation x p n (x) = a n p n+1 (x) + b n p n (x) + c n p n−1 (x), vanishes. In this case, discrete-continuous Sobolev orthogonality must be considered since, due to Favard's Theorem (see [6], Theorem 4.4), there is no regular linear functional associated with these families of polynomials.
The main strategy for dealing with Jacobi and Laguerre polynomials with non-standard parameters (Gegenbauer polynomials included) is to use the fact that the derivatives of the classical orthogonal polynomials are again polynomials of the same family with a shift in their parameters. For instance, in [15], the Sobolev orthogonality for the Laguerre polynomials  n (x)} n≥0 and the rest of the functionals being derivatives of Dirac deltas supported on the origin of the real line. Later, in [18], a unified approach to the orthogonality of the generalized Laguerre polynomials {L (α) n (x)} n≥0 was given for any real value of the parameter α, by proving their orthogonality with respect to a more general Sobolev non-diagonal inner product that reduces to (1.1) when −α = N ∈ N.
In [3], a general approach was given for dealing with the orthogonality of Gegenbauer polynomials {C (−m+1/2) n (x)} n≥0 with non standard parameters, that is, for m ∈ N. Therein, the bilinear form (1.1) is expressed as the sum of two bilinear forms: a discrete form D (u k , 0 ≤ k ≤ 2 m, are derivatives of Dirac deltas supported on the endpoints of the interval of orthogonality [−1, 1]), and a continuous form C defined by means of C( p, q) = u 2m , p (2m) q (2m) , where u 2m is the regular linear functional associated with the Gegenbauer polynomials {C (m+1/2) n (x)} n≥0 . We must remark that the strategy followed in [3] is to prove the existence of a 2m × 2m symmetric positive definite matrix A such that D( p, q) = P A Q t , (1.3) where . . , f (m−1) (−1)), and, therefore, B( p, q) = P A Q t + u 2m , p (2m) q (2m) .
Observe that B acts on polynomials of degree ≤ 2m −1 only through D, and, for polynomials of degree ≥ 2m, B acts only through C. In the same way, in [19], the Jacobi polynomials {P were proven to be orthogonal with respect to . . , f (N −1) (1)), A is a symmetric positive definite matrix and u N is the regular linear functional associated with the Jacobi polynomials {P (0,β+N ) n (x)} n≥0 . The study of the Jacobi polynomials {P (−k,− ) n (x)} n≥0 , k, ∈ N, presents an additional challenge. As shown in [20], the hypergeometric expression for P (α,β) n (x) is valid for arbitrary complex values of the parameters α and β. Nevertheless, a reduction of the degree of P (α,β) n (x), n ≥ 1, occurs if and only if n = −α − β − j for a certain integer 1 ≤ j ≤ n. In [1], this issue is dealt with by defining the monic generalized Jacobi polynomial P (−k,− ) n (x) by means of the following formula: where h = k + −n −1 if (k + )/2 ≤ n < max{k, }, and h = n otherwise. The polynomials defined in (1.4) inherit several of the important properties of the Jacobi polynomials and satisfy deg P (−k,− ) n = n, n ≥ 0. Moreover, these polynomials are orthogonal with respect to the bilinear where C is a discrete bilinear form as in (1.3) with A being an N × N symmetric matrix, and u is the regular linear functional associated with {P ( ,k) n (x)} n≥0 . Using some ideas from [17], a construction of discrete-continuous bilinear forms is presented in [19], that is valid for general orthogonal polynomials { p n (x)} n≥0 , not necessarily the Jacobi and Laguerre polynomials. Therein, the discrete part D is defined by D( p n , p m ) = k n δ n,m , n, m ≤ N − 1, with arbitrary positive constants k n . Therefore, D acts on the linear space of polynomials of degree at most N − 1, and its Gram matrix associated with the basis { p n (x)} N −1 n=0 is a diagonal matrix K = diag{k 0 , k 1 , . . . , k N −1 }. A non singular matrix H is defined as the inverse of a matrix whose entries are the polynomials p 0 (x), p 1 (x), . . . , p N −1 (x) and possibly their derivatives evaluated at the roots of polynomials p N (x) (taking into account their multiplicities) which are known. Then, where P and Q are vectors whose entries are the evaluations of p(x) and q(x), and possibly their derivatives, at the roots of p N (x). Note that A = H K H t is the Gram matrix of D with respect to the basis of Lagrange interpolation polynomials associated with the roots of p N (x). In this way, D is the restriction of the bilinear form B to the linear space of polynomials of degree at most N −1, and, thus, B = D +C for some appropriate continuous bilinear form C.
Finally, [16] studies the orthogonality of the polynomial solutions of the linear secondorder differential equation where φ and ψ are non zero polynomials of degree at most 2, and λ n are constants. Therein, we find the following result.
The differential operator L is symmetric with respect to B, that is, for all polynomials p, q. (c) The linear functionals u 0 and u 1 satisfy the distributional differential equations Using the fact that the classical orthogonal polynomials are characterized by satisfying the differential equation (1.5) (see [4] as well as [16,Theorem 3.5.]), [16] deals with some examples of classical orthogonal polynomials with non standard parameters, including the Laguerre polynomials with parameter α = −1; and the Jacobi polynomials with parameters α = β = −1, α = −1 and −β / ∈ N, and −α / ∈ N and β = −1. The general aim of this paper is motivated by Theorem 3.3 of [16] in the sense that we study the orthogonality of polynomial solutions of (1.5) when L is symmetric with respect to (1.1) for arbitrary N ≥ 0. In this way, we show that the discrete-continuous Sobolev orthogonality of classical orthogonal polynomials with non standard parameters arises from the symmetry of the operator L and, in fact, we recover the matrix representation of the discrete part of the corresponding bilinear forms. In contrast with [1], we deal with the reduction in degree of the Jacobi polynomials P (−k,− ) n (x), k, ∈ N, by explicitly constructing a polynomial solution of (1.5) of degree n for each n ≥ 0.
In Sect. 2, we present basic facts about orthogonal polynomials and, in particular, the classical Hermite, Laguerre, Bessel, and Jacobi polynomials. The main results of this paper are collected in Sect. 3, where we study conditions for a linear second-order differential operator to be symmetric with respect to the bilinear form (1.1). We show that the linear functionals in B must satisfy the related distributional Pearson equations. Hence, in this section we also study the solutions of the Pearson equations. Lastly, in Sect. 4, we apply our results to deduce the Sobolev orthogonality for polynomial solutions of the differential equations corresponding to the classical orthogonal polynomials with non-standard parameters.

Preliminaries
In this section, we present the basic facts needed to establish our main results.

Orthogonal polynomials
We denote by R the set of real numbers, and N denotes the set of positive integers. For each integer n ≥ 0, we denote by Π n the linear space of real polynomials of degree at most n, and let Π = n≥0 Π n be the collection of all such polynomials.
We call any linear functional u : Π → R a moment functional, and the image of a polynomial p ∈ Π under u will be denoted by u, p . We denote by Π the linear space of moment functionals, that is, Π = {u : Π → R | u is linear}. Observe that Π is the algebraic dual of Π. Given a sequence of real numbers {μ n } n≥0 , a moment functional u can be defined by means of its moments as u, x n = μ n , n ≥ 0, and extended by linearity to all of Π.
Let D : Π → Π denote the distributional differential operator. Given a moment functional u, its derivative Du is the moment functional defined by Du, p = − u, p , ∀ p ∈ Π, and the left multiplication of u by a polynomial q ∈ Π is the moment functional qu defined by qu, p = u, q p , ∀ p ∈ Π.
We call any map B : Π × Π → R a bilinear form if for all p, q, r ∈ Π and λ ∈ R, it holds that B(λ p + q, r ) = λ B( p, r ) + B(q, r ) and B(r , λ p + q) = λ B(r , p) + B(r , q), that is, B is linear in each one of its entries. Additionally, if B( p, q) = B(q, p) for all p, q ∈ Π, then we say that B is symmetric.
Let { p n (x)} n≥0 be a sequence in Π such that deg p n = n for n ≥ 0. Then { p n (x)} n≥0 is a basis for Π. In addition, if B( p n , p m ) = 0 for n = m and B( p n , p n ) = 0 for n ≥ 0, we say that { p n (x)} n≥0 is an orthogonal polynomial sequence (OPS) associated with B.
Given a symmetric bilinear form B, there is not always an OPS associated with B. If an OPS associated with B exists, then B is called quasi-definite. A bilinear form B is positive definite if B( p, p) > 0 for every non-zero polynomial p ∈ Π. Positive definite bilinear forms are quasi-definite and, thus, an OPS associated with B exists.
Following [16], the quasi-definite and positive definite character of a symmetric bilinear form B can be characterized in terms of the principal minors of the Gram matrix defined from the elements μ i, j : (resp. Δ n (B) > 0, for all n ≥ 0). Observe that, since B is symmetric, μ i, j = μ j,i , i, j ≥ 0, and, thus, Δ n (B), n ≥ 0, are the principal minors of a symmetric Gram matrix. Given u ∈ Π , we can define the bilinear form B( p, q) = u, p q , p, q ∈ Π. Then, we will say that u is quasi-definite (resp. positive definite) if and only if B is quasi-definite (resp. positive definite), and that there is an OPS associated with u. In this case, B(x i , x j ) = u, x i+ j and, thus, Δ n (B), n ≥ 0, are the principal minors of a Gram matrix with Hankel structure.

Classical moment functionals
A quasi-definite moment functional u is called classical if there are non-zero polynomials φ and ψ, with deg φ ≤ 2 and deg ψ = 1, such that u satisfies the distributional Pearson equation or, equivalently, An OPS associated with a classical moment functional is called a classical OPS.
Although there are several properties that characterize the classical moment functionals, we focus our attention on two of them: the first one given in 1929 by Bochner [4] and the second one given in 1935 by Hahn [11]. We state both characterizations in the following theorem.
Theorem 1 Let u be a quasi-definite functional and {P n (x)} n≥0 an OPS associated with u. The following statements are equivalent.

[11]
There is a non zero polynomial φ with deg φ ≤ 2, such that P n+1 n≥0 is a sequence of orthogonal polynomials associated with the moment functional v = φ u.
It is well known (see [4] as well as [12,17]) that, up to affine transformations of the independent variable, the only families of classical orthogonal polynomials are the Hermite, Laguerre, Jacobi, and Bessel polynomials. The classical orthogonal polynomials are classified according to the canonical form of the polynomial φ(x).
We remark that when φ(x) ≡ 0, the polynomial solutions of (2.1) are the elements of the sequence {x n } n≥0 , which can not be associated with a quasi-definite moment functional.
The following explicit expressions for the Hermite, Laguerre, and Jacobi polynomials appear in [20], and the explicit expression for the Bessel polynomials is found in [13, p. 108]. Hermite: where, as usual, denotes the Pochhammer symbol, and denotes the extended binomial coefficient. As shown in [20], the expression for P (α,β) n (x) is valid for arbitrary complex values of the parameters α and β. Nevertheless, a reduction of the degree of P n (x) is valid for arbitrary values of the parameter α (see [20]), and no reduction in the degree ever occurs. For −a ∈ N, the expression for B For arbitrary values of the parameters, the above expressions for H n , L n are polynomials solutions of the differential equation (2.1) with the corresponding polynomials φ(x) and ψ(x). Namely, they satisfy the following linear second-order differential equations: n : It is important to note that if either α, β, or α + β + 1 are negative integers, then the Jacobi polynomials {P are associated with a quasi-definite moment functional u α,β satisfying the Pearson equation Furthermore, for α, β > −1, u α,β is positive definite and has the integral representation n (x)} n≥0 are associated with a quasi-definite moment functional u α satisfying the Pearson equation The polynomials {L (−m) n (x)} n≥0 , with m ∈ N can not be associated with a quasi-definite moment functional.
The Hermite polynomials {H n (x)} n≥0 are associated with the positive definite moment functional u defined by and satisfying the Pearson equation n (x)} n≥0 are only associated with a non-positive definite moment functional u (a) defined by where c is the unit circle oriented in the counter-clockwise direction, satisfying the Pearson equation

Sobolev bilinear forms and differential operators
From Theorem 1, we have that a family of classical orthogonal polynomials associated with a quasi-definite moment functional u are solutions of the linear second-order differential equation (2.1) and are simultaneously orthogonal with respect to the symmetric Sobolev bilinear form where φ is a non zero polynomial of degree at most 2. In this section, we study the orthogonality of polynomial solutions of (2.1) with respect to a Sobolev bilinear form involving derivatives up to a fixed but arbitrary order. Fix an integer N ≥ 0, and let u k , 0 ≤ k ≤ N be moment functionals in Π . A symmetric Sobolev bilinear form can be defined in the following way: (3.1) We say that two polynomials p, q ∈ Π are orthogonal with respect to B N if We also define the linear second-order differential operator L acting on Π as where φ and ψ are polynomials such that deg φ ≤ 2 and deg ψ ≤ 1. Its formal Lagrange adjoint L * is the operator acting on Π satisfying and, thus, In the sequel, we will be interested in the polynomial solutions { p n (x)} n≥0 of the differential equation Remark 1 It was shown in [16] that (3.3) is admissible if and only if it has a unique linearly independent polynomial solution of degree n for all n ≥ 0, and It is straightforward to verify that the differential equation is admissible if and only if, for all n ≥ 0, Definition 1 We say that the operator L is symmetric with respect to B N if The symmetry of the operator L and an additional hypothesis imply that the polynomial solutions of the differential equation (3.3) are orthogonal with respect to B N . Theorem 2 Let L be the differential operator defined in (3.2), and assume that it is symmetric with respect to the Sobolev bilinear form (3.1). If p, q ∈ Π are polynomials satisfying that is, p and q are orthogonal with respect to B N .
Proof From the symmetry of L with respect to B N , we have Corollary 1 Let L be the differential operator defined in (3.2). Let { p n (x)} n≥0 be a sequence of polynomials satisfying the differential equation and assume that the differential equation is admissible.
If L is symmetric with respect to B N , then The following result provides a partial reciprocal of the previous corollary. (3.1). Assume that B N is quasi-definite, and denote by { p n (x)} n≥0 an OPS associated with B N . Let L be the operator defined in (3.2), and assume that d n := 1 2 n φ + ψ = 0, for all n ≥ 0. Then, the following statements are equivalent.

Theorem 3 Let B N be the bilinear form defined in
(i) The operator L is symmetric with respect to B N .
(ii) For each n ≥ 0, there is a constant λ n such that L[ p n ] = λ n p n .
Proof (i)⇒(ii) For n ≥ 0, observe that since d n = 0, the operator L preserves the degree, thus L[ p n ] is a polynomial of degree n. Therefore, Using the symmetry of L, we get Then, L[ p n ] = c n p n and (ii) holds with λ n = c n .
(ii)⇒(i). Since B N is bilinear, it suffices to show (i) for a basis of polynomials. For n, m ≥ 0, we have Clearly, (i) follows.
We are interested in deducing the necessary and sufficient conditions for the differential operator (3.2) to be symmetric with respect to B N in (3.1). With this in mind, we state the following preliminary result. Lemma 2 Let n ≥ 0, and let L be the operator defined in (3.2). Then, for all p ∈ Π and all u ∈ Π , the following relation holds, Proof For n = 0, it can be easily verified that, for all p ∈ Π and all u ∈ Π , For n ≥ 1, using Leibniz product rule and the fact that deg φ ≤ 2 and deg ψ ≤ 1, we compute Similarly, Using (3.4) and Now, we are ready to state the main result of this section. (3.1), and let L be the second-order differential operator defined in (3.2). Then L is symmetric with respect to B N if and only if the moment functionals u k , 0 ≤ k ≤ N , satisfy the Pearson equations

Theorem 4 Let B N be the bilinear form defined in
Proof Clearly, L is symmetric with respect to B N if and only if, for all p, q ∈ Π, Then, L is symmetric with respect to B N if and only if, for all p ∈ Π, (3.6) On one hand, assume that, for 0 ≤ k ≤ N , u k satisfies D(φ u k ) = (ψ + k φ ) u k . Then, (3.6) holds for all p ∈ Π and therefore L is symmetric with respect to B N .
On the other hand, assume that (3.6) holds for all p ∈ Π. Then, letting p(x) = 1, (3.6) reads From Lemma 1, we deduce that D(φ u 0 )−ψ u 0 = 0. Consequently, the term corresponding to k = 0 in (3.6) is 0, and we can sum from k = 1 to k = N . Now, letting p(x) = x, (3.6) reads Using Lemma 1 twice, we deduce that D(φ u 1 ) − (ψ + φ ) u 1 = 0, and, consequently, the summation in (3.6) goes from k = 2 to k = N . Repeating this process for 3 ≤ m ≤ N by We turn our attention to the solutions of the Pearson equations in Theorem 4. Writing and substituting in the previous equation, we obtain and (3.7) follows. It is easy to verify that the implication in the opposite direction holds by inverting each of the previous steps. (ii) If d n = 0, n ≥ 0, then (3.7) is a linear recurrence relation. It follows that (3.7) is uniquely solvable for μ n+1 , n ≥ 0, once μ 0 is fixed. Then, the moment functional u satisfying D(φ u) = ψ u is uniquely determined up to a constant factor.

Remark 2
Observe that if the differential equation

Sobolev orthogonal polynomial solutions of differential equations
In this section, we focus our attention on the differential equations satisfied by the Hermite, Laguerre, Jacobi, and Bessel polynomials. In particular, we study in great detail the Sobolev orthogonality of the polynomial solutions, including the case when the parameters in the Jacobi and Laguerre polynomials have negative integer values.
As we have already mentioned in Sect. 2, though valid for arbitrary complex values of the parameters, the expression for P = n for all n ≥ 0. Nevertheless, it is shown in the sequel that there is a polynomial sequence { p n (x)} n≥0 with deg p n = n satisfying this differential equation.
The following proposition will be used to construct polynomial solutions of this differential equation and, in general, of (3.3).

Proposition 3
Let L be the differential operator defined in (3.2). Assume there is a non zero polynomial ρ satisfying the differential equations

2)
where h is a polynomial of degree ≤ 1. Let k = deg ρ. Moreover, let { p n (x)} n≥0 be a sequence of polynomials satisfying the differential equation Then the polynomials q n (x) = ρ(x) p n−k (x), n ≥ k, satisfy the differential equation Proof Using (4.1) and (4.2), we obtain

Laguerre polynomials
As it was explained in Sect. 2, for arbitrary values of the parameter α, the Laguerre polynomials {L Define the linear second-order differential operator Let N be a fixed non-negative integer. By Theorem 4, L (α) is symmetric with respect to the Sobolev bilinear form B N defined in (3.1) if and only if the moment functionals u k , 0 ≤ k ≤ N , satisfy the Pearson equations In this case, (3.8) reads d n = −1, n ≥ 0.
We will proceed to study the Sobolev orthogonality for Laguerre polynomials with standard and non-standard parameters, that is, for −α / ∈ N and −α ∈ N, respectively.

Sobolev orthogonality for standard parameter
For −α / ∈ N, u 0 can be chosen to be the quasi-definite moment functional defined as Moreover, from Proposition 2, we deduce that u k = a k x k u 0 , 1 ≤ k ≤ N , where a k are constants. Hence, (3.1) reads where a 0 = 1. Since L (α) is symmetric with respect to B N , the Laguerre polynomials satisfy   Therefore, we get

Sobolev orthogonality for non standard parameter
where b j,ν are constants. Using (4.3), we get where the constants b j = b j,0 are free parameters. Then (3.1) reads If we define the discrete bilinear form D as and the continuous bilinear form C as

then (4.4) can be expressed as
Additionally, observe that the polynomial ρ(x) = x m satisfies the differential equations

Remark 3
The matrix expression of the discrete bilinear form D can be written as follows. For 0 ≤ j ≤ m − 1, let F( j) be the symmetric matrix of size (m − j) × (m − j) with entries given by and, for 0 ≤ k ≤ m − 1, let G(k) be the symmetric matrix of size m × m defined as G(0) = F(0), and Then, Moreover, let L(0) be the non-singular lower triangular matrix of derivatives of Laguerre polynomials given by Then, we have the m × m diagonal matrix Clearly, if B N is quasi-definite, then H must be chosen to be a positive definite matrix.

Jacobi polynomials
For arbitrary complex values of α and β, the Jacobi polynomials {P (α,β) n (x)} n≥0 are polynomial solutions of the differential equation x] y = λ n y, λ n = −n(n + α + β + 1). (4.5) Define the differential operator and fix a non negative integer N . We know from Theorem 4 that L (α,β) is symmetric with respect to the Sobolev bilinear form B N defined in (3.1) if and only if the moment functionals u k , 0 ≤ k ≤ N , satisfy the Pearson equations In this case, (3.8) reads Here, we study the Sobolev orthogonality for the polynomial solutions of (4.5) when the values of the parameters fall in one of the following four cases:

Sobolev orthogonality for standard parameters
For −α, −β, −(α +β +1) / ∈ N, let u 0 ∈ Π be the quasi-definite moment functional defined as and, for 1 ≤ k ≤ N , let u k = a k (1 − x 2 ) k u 0 , where a k are constants. Then, B N in (3.1) with the moment functionals u k , 0 ≤ k ≤ N , reads with a 0 = 1. Moreover, since each moment functional u k satisfies the Pearson equation (4.6), the operator L (α,β) is symmetric with respect to B N .
From the symmetry of L (α,β) with respect to B N , we have that

Sobolev orthogonality for type I non-standard parameters
For −α = m ∈ N and −β / ∈ N, we have that d n+2 m = −(n + m + β + 2) = 0, n ≥ 0. and that the differential equation satisfied by the Jacobi polynomials {P is uniquely solvable up to a constant factor. Then u m can be chosen to be the quasi-definite moment functional associated with the Jacobi polynomials {P (0,β+m) n (x)} n≥0 , that is, By Proposition 2, where a i are constants (a 0 = 1).
For 1 ≤ j ≤ m, let the moment functional u m− j satisfy the Pearson equation or, equivalently, Hence, we deduce that where c j is a free parameter. Taking u k , 0 ≤ k ≤ N , as above, (3.1) reads with a 0 = 1. Since the Jacobi polynomials {P (−m,β) n (x)} n≥0 satisfy the admissible differential equation (4.8), then it follows from the symmetry of L (−m,β) with respect to B N and Corollary 1 that Observe that B N can be expressed as the sum of two bilinear forms: a discrete bilinear form and a continuous bilinear form In this way

Remark 4
The discrete bilinear form D I can be written in matrix form. Indeed, for 1 ≤ j ≤ m, let A( j) be the symmetric matrix of size j × j with entries given by where 0 is the zero matrix of appropriate size. Then, Moreover, taking into account that D I coincides with the restriction of B N to the linear space of polynomials of degree at most m − 1, we deduce that and P(1) is the non-singular matrix of derivatives of Jacobi polynomials (1) n,k=0,1,...,m−1 , given by Observe that the numbers h (−m,β) n must depend on the free parameters c j . Therefore, given an arbitrary m × m diagonal matrix H, D I can be explicitly constructed by means of Clearly, if B N is quasi-definite, then H must be chosen to be a positive definite matrix.

Sobolev orthogonality for type II non-standard parameters
The Jacobi polynomials satisfy the important identity ([20, p. 59]) Then it is possible to give the Sobolev orthogonality for the polynomial sequence {P Observe that for −α / ∈ N and −β = ∈ N, the differential equation is admissible and d n+2 = −(n + α + + 2) = 0, n ≥ 0.
Fix N ≥ . Constructing the solutions of the Pearson equations (4.6) as in the previous subsection, we obtain with b i being a constant (b 0 = 1), and for 1 ≤ j ≤ , where e j is a free parameter. Therefore, (3.1) with the above moment functionals is written as where b 0 = 1. Since the Jacobi polynomials {P (α,− ) n (x)} n≥0 satisfy (4.12), then it follows from the symmetry of L (α,− ) with respect to B N and Corollary 1, that they satisfy If we define the discrete bilinear form D I I by and the continuous bilinear form C I I by Moreover, using (4.9) and (4.11), we obtain Remark 5 Analogously to Remark 4, for 1 ≤ j ≤ , let C( j) be the symmetric matrix of size j × j with entries given by and, for 0 ≤ k ≤ −1, let D(k) be the symmetric matrix of size × defined as D(0) = C( ), and ), 0 ≤ n ≤ − 1, and P(−1) is the nonsingular matrix of derivatives of the Jacobi polynomials given by Moreover, since the numbers h (α,− ) n depend on the free parameters e j , D I I can be explicitly constructed by means of with H being an arbitrary × diagonal matrix. If B N is quasi-definite, then H must be chosen to be a positive definite matrix.

Sobolev orthogonality type III non-standard parameters
Here, we consider the polynomial solutions of (4.5) with −α = m and −β = , m, ∈ N. In this case, the differential equation is not admissible. Moreover, the polynomial P (−m,− ) n (x) suffers a reduction in degree whenever n = m + − k for some integer 1 ≤ k ≤ n. Hence, we are concerned with constructing polynomial solutions {q n (x)} n≥0 of (4.5) such that deg q n = n, and also to give the Sobolev orthogonality for these solutions.
First, we will construct polynomial solutions of (4.5). In [1], the Jacobi polynomials were redefined in order to circumvent the reduction in degree of P (−m,− ) n (x), and there it was shown that these redefined Jacobi polynomials satisfy (4.5) and are orthogonal with respect to a Sobolev bilinear form. We will not use this alternative definition of the Jacobi polynomials, but instead we use Theorem 3 to construct polynomial solutions of (4.5).
Due to (4.11), we can assume without any loss of generality that ≤ m. Observe that for 0 ≤ n ≤ − 1, the polynomial P (−m,− ) n (x) has degree n and satisfies (4.5). No degree reduction occurs in this case.
For ≤ n ≤ m + − 1, the polynomial q n (x) = (1 + x) P (−m, ) n− has degree n and satisfies (4.5). Indeed, let ρ 1 (x) = (1 + x) . The polynomial ρ 1 satisfies the differential equations Moreover, the polynomial P For n ≥ m + , define the polynomialsq n (x) = (1 − x) m (1 + x) P (m, ) n−m− (x). Clearly, degq n = n. We will show that these polynomials satisfy (4.5). Define the polynomial (1 + x) , and observe that ρ 2 satisfies the differential equations Taking into account that the Jacobi polynomials P Altogether, the sequence of polynomials {P Observe that d n+2 (m+ ) = 0 for n ≥ 0. By Proposition 2, we have where a i are constants (a 0 = 1), and u m+ satisfies Then u m+ can be chosen to be the quasi-definite moment functional associated with the Jacobi polynomials {P ( ,m) n (x)} n≥0 , that is, For 1 ≤ j ≤ m + , let the moment functional u m+ − j satisfy the Pearson equation For 1 ≤ j ≤ m, we obtain where c j is a free parameter. Similarly, for 1 ≤ r = j − m ≤ , we get where e r is a free parameter. Taking the moment functionals u k , 0 ≤ k ≤ N , as above, the bilinear form (3.1) reads (4.15) Let D I I I be the discrete bilinear form defined by and let C I I I be the continuous bilinear form defined by Then (4.15) can be written as On the other hand, we compute Remark 6 Using (4.10) and (4.13), the bilinear form D I I I can be rewritten as If the Sobolev bilinear form B N is quasi-definite, then H must be chosen to be a positive definite matrix.

Sobolev orthogonality for type IV non-standard parameters
Now, the case when −α = −β = m ∈ N is considered. As we will see here, this is not a particular case of the third type of non-standard parameters studied above, in the sense that the continuous part of the Sobolev bilinear form is not directly deduced by setting m = in Hence, in general, C I I I and C I V do not coincide. Now, we proceed to find solutions to the remaining Pearson equations. For 1 ≤ j ≤ m, let the moment functional u m− j satisfy the Pearson equation There are two linearly independent solutions. Namely, u m− j,1 and u m− j,2 given by Then, where c j and e j are free parameters. Therefore, if we define the discrete bilinear form D I V as D I I I with m = , that is, Observe that the polynomial sequence (4.14) with m = satisfies the differential equation

Hermite polynomials
The Hermite polynomials {H n (x)} n≥0 constitute a sequence of polynomial solutions of the the admissible equation L[y] ≡ y − 2 x y = −2 n y.
By Theorem 4, we have that for a non-negative integer N , the differential L is symmetric with respect to the bilinear form B N defined in (3.1) if and only if D(u k ) = −2x u k for 0 ≤ k ≤ N . Thus, if we take u 0 such that then from Proposition 2, u k = a k u 0 and therefore (3.1) reads where a k are constants (a 0 = 1). From the symmetry of L with respect to B N and Corollary 1, we have B N (H n , H m ) = 0 for n = m.

Bessel polynomials
Let L (a) be the differential operator defined by Recall from Sect. 2 that for arbitrary complex values of a, the Bessel polynomials {B In this case, (3.8) reads d n = n + a, n ≥ 0. For −a / ∈ N, we can choose u k = a k x 2 k u 0 , 0 ≤ k ≤ N , where u 0 , p = 1 2π i c p(z) z a−2 e −2/z dz, ∀ p ∈ Π, c is the unit circle oriented in the counter-clockwise direction, and a k are constants (a 0 = 1). Then, we have that L (a) is symmetric with respect to the bilinear form In this case, deg B (a) n = n, n ≥ 0. Moreover, the differential equation satisfied by the Bessel polynomials is admissible. We know from Corollary 1 that where ν is the smallest integer bigger that ν. Unfortunately, in this case the differential equation (4.17) does not have a polynomial solution of degree n for m 2 + 1 ≤ n ≤ m + 1. Hence, we do not give a Sobolev orthogonality for this case.