Recurrence relations for orthogonal polynomials for PDEs in polar and cylindrical geometries

This paper introduces two families of orthogonal polynomials on the interval (−1,1), with weight function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega (x)\equiv 1$$\end{document}ω(x)≡1. The first family satisfies the boundary condition \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p(1)=0$$\end{document}p(1)=0, and the second one satisfies the boundary conditions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p(-1)=p(1)=0$$\end{document}p(-1)=p(1)=0. These boundary conditions arise naturally from PDEs defined on a disk with Dirichlet boundary conditions and the requirement of regularity in Cartesian coordinates. The families of orthogonal polynomials are obtained by orthogonalizing short linear combinations of Legendre polynomials that satisfy the same boundary conditions. Then, the three-term recurrence relations are derived. Finally, it is shown that from these recurrence relations, one can efficiently compute the corresponding recurrences for generalized Jacobi polynomials that satisfy the same boundary conditions.

order. Generalization to other (non-integer) indices was carried out in Guo et al. (2009) to obtain families of orthogonal polynomials for Chebyshev spectral methods or problems with singular coefficients. However, although these GJPs can be described in terms of short linear combinations of Legendre polynomials, at least for certain index pairs of interest (Guo et al. 2009;Shen 2003), the three-term recurrence relations characteristic of families of orthogonal polynomials have not been developed in these cases.
In this paper, we use the bases from Shen (1997) to develop families of polynomials that are orthogonal with respect to ω(x) ≡ 1 and satisfy the requisite boundary conditions, to facilitate transformation between physical and frequency space without using functions such as the Legendre polynomials that lie outside of the solution space. These families can also be efficiently modified to work with alternative weight functions, thus leading to the development of new numerical methods. In particular, it is demonstrated that these new families can be used to obtain three-term recurrence relations for the GJPs that satisfy the same boundary conditions. The outline of the paper is as follows. In section "Variational formulation", we provide context for these families of polynomials by adapting the variational formulation employed in Shen (1997) to the time-dependent PDE (1)-(3). In section "The case m = 0" we develop orthogonal polynomials with unit weight function satisfying the boundary conditions p(1) = 0. In section "The case m ≠ 0" we do the same for the boundary conditions p(−1) = p(1) = 0. In section "Recurrence relations for generalized jacobi polynomials" we describe how these families of orthogonal polynomials can be efficiently modified to obtain three-term recurrence relations for GJPs as described in Guo et al. (2009), Shen (2003. Concluding remarks and directions for future work are given in section "Conclusions".

Variational formulation
In this section, we describe one possible context in which the sequences of orthogonal polynomials discussed in this paper can be applied.

Conversion to polar coordinates
We consider the reaction-diffusion equation on a unit disk where α is a constant.
Following the approach used in Shen (1997) for a steady-state problem, we can convert the IBVP in (1)-(3) to polar coordinates by applying the polar transformation x = r cos θ , y = r sin θ and letting u(r, θ ) = U (r cos θ , r sin θ ), f (r, θ ) = F (r cos θ, r sin θ ). The resulting problem in polar coordinates is as follows: The solution is represented using the Fourier series The Fourier coefficients u 1,m (r, t), u 2,m (r, t) must satisfy the boundary conditions u 1,m (1, t) = u 2,m (1, t) = 0 for m = 0, 1, 2, . . . Due to the singularity at the pole r = 0, we must impose additional pole conditions on (5) to have regularity in Cartesian coordinates. For u(r, θ , t) to be infinitely differentiable in the Cartesian plane, the additional pole conditions are Shen (1997) By substituting the series (5) into (4) and applying the pole conditions in (6), we obtain the following ODEs, for each nonnegative integer m: where u and f are now generic functions.

Legendre-Galerkin method
To approximate (9) using the Legendre-Galerkin method, we let ω = 1 and we have to where I N is the interpolation operator based on the Legendre-Gauss-Lobatto points. That is, and L N is the Legendre polynomial of degree N.

The case m = 0
In the case where m = 0, (10) reduces to As before, we let L k (t) be the kth-degree Legendre polynomial, and define X N (0) to be the space of all polynomials of degree less than or equal to N that vanish at 1. This space can be described as Shen (1997) where φ i (t) is the ith basis function. By applying the Gram-Schmidt process (Burden, Faires 2005) to these basis functions, φ i (t), we can obtain a new set of orthogonal polynomials that will be denoted by φ i , i = 0, 1, 2, . . ., where the degree of φ i and φ i is i + 1. The new basis functions, φ i , can be found by computing due to the orthogonality of the Legendre polynomials, thus greatly simplifying the computation of φ i . (10) To start the sequence {φ i }, we let so then and The first several polynomials φ 0 ,φ 1 , . . . ,φ 4 are shown in Fig. 1. Now, comparing φ 1 with φ 1 and φ 2 with φ 2 , we can find a general formula for the φ i in terms of φ i . By subtracting φ i from φ i , we obtaiñ x 2 + 11 6 x − 1 6 . and This suggests a simple recurrence relation for φ i in terms of φ i . Before we prove that this relation holds in general, we need the following result.

Lemma 1 Let
Proof We proceed by induction. For the base case, we have For the induction step, we assume that there is a k > 0, such that N k−1 = 2(k+1) 2 k 2 (2k+1) . We must show that the formula found in Eq. (15) is true for k. Given φ k = φ k + k k+1 2φ k−1 , and using (13) We can now establish the pattern seen in (13), (14).
Proof Again we proceed by induction. For the base case, we will show that the theorem holds when i = 1: Note that Eq. (18) is equivalent to Eq. (13). For the induction step, we assume that there is a j ≥ 0, such that Richardson andLambers SpringerPlus (2016) 5:1567 We show that (17) holds when i = j + 1. We have Therefore, using Lemma 1 and (16), we obtain We now prove a converse of Theorem 1.
All orthogonal polynomials satisfy a general three-term recurrence relation that has the form where α j , β j and γ j are constants. By enforcing orthogonality, we obtain the formulas First, we will find the value of α j .
Theorem 3 Let α j be defined as in (21).
Proof Base case: When j = 0, we use (21) to obtain For the induction hypothesis, we assume there is a j > 0 such that , we obtain Now, from the recurrence relation for Legendre polynomials, we obtain .

and
To calculate the middle term in Eq. (24) we will multiply 2c j by the result from Eq. (26): We rearrange the formula for α j−1 to obtain the following: Therefore, 2c j φ j−1 , xφ j = 2 j j + 1 2 2 j 2 + 3j + 3 Now we can use the results from Eqs. (25)-(28) to determine the numerator of α j .
Hence, Now, we will find the value of β k .
Proof For the base case, we consider j = 0: For the induction step, we assume there is a j ≥ 0 such that β j−1 = j+1 2j+1 . From Using the recurrence relation for Legendre polynomials, we obtain and We then have The last term in (29) is obtained as follows:

We then have
We rearrange the formula for β j−1 to obtain the following: Therefore, Now we can use the results from Eqs. (30)-(33) to determine the numerator of β j .
In summary, the polynomials {φ i } satisfy the recurrence relation We can rewrite Eq. (19) as φ j − c jφj−1 = φ j . In matrix form, we have , with x being a vector of at least n + 2 Legendre-Gauss-Lobatto points. This ensures that the columns of Φ are orthogonal.
Then, given f ∈ X n+1 (0), we can obtain the coefficients f i in by simply computing f i = �φ i , f �/N i , where N i is as defined in (15). Then the coefficients f i in can be obtained by solving the system Cf =f using back substitution, where C is as defined in (36). These coefficients can be used in conjunction with the discretization used in Shen (1997), which makes use of the basis {φ i }.

The case m � = 0
In the case where m � = 0, we work with the space As discussed in Shen (1997), this space can easily be described in terms of Legendre polynomials: Applying the Gram-Schmidt process to the basis functions {φ i }, we obtain a new set of orthogonal polynomials that will be denoted as {φ i }. These basis functions are obtained in the same way as in Eq. (11). First, we let and Then, we have and The graphs of the first several members of the sequence {φ i } are shown in Fig. 2. Again, we will compare φ 2 with φ 2 and φ 3 with φ 3 to find a general formula for the values of φ i . We obtain the following formulâ x. (37) and These results suggest a simple recurrence relation for φ i in terms of φ i and φ i−2 , in which the coefficient of φ i−2 is a ratio of triangular numbers d i = i(i − 1)/[(i + 1)(i + 2)]. We therefore define with initial conditions To prove that these polynomials are actually orthogonal, we first need this result.

∀j ≥ 2 and
Proof For the base case, we compute N 0 and N 1 directly. We have and For the induction step, we assume there is a j > 2 such that N j−2 = 2j(j+3) j(j+1)(2j+3) . Now, we must show that the formula (40) is true for j. We have Theorem 5 Let φ i be obtained by orthogonalizing φ i against φ 0 ,φ 1 , . . . Then φ 0 = φ 0 , φ 1 = φ 1 , and Proof For the base case, we first show that φ 1 = φ 1 and φ 0 = φ 0 are already orthogonal. We have Next, we show directly that the theorem holds when j = 2: For the induction step, we assume that φ 0 , . . . ,φ j−1 are all orthogonal, where j ≥ 2, and that Using Lemma 2, we obtain We now confirm that the polynomials defined using the recurrence (41) are orthogonal.
Case 3: j = k − 1. If k ≥ 3, then we have If k = 2, then the steps are the same, except that the term with φ k−3 is not present.
Like all families of orthogonal polynomials, the {φ k } satisfy the recurrence relation By analogy with (21), (22) and (23), we have Because φ j contains only terms of odd degree if j is odd and of even degree if j is even, just like the Legendre polynomials, it is easily shown that α j = 0 for j = 1, 2, . . . We will now find the values of β j and γ j .
Proof We show the base case j = 0 directly: For the induction step, we assume there is a j ≥ 0 such that β j−1 = j+2 2j+3 .
Then, using (45), we have β j = φ j+1 ,xφ j φ j+1 ,φ j+1 . For the numerator, we have We now compute each part of this numerator as follows: Then For the third term in (47), we have j + 1 2j + 1 2j + 5 , and therefore We rearrange the formula for β j−2 to obtain the following: Therefore, Now we can use the results from Eqs. (48)-(51) to determine the numerator of β j .
we finally obtain and delete the last row and column to obtain Ĵ n−2 .
To overcome the obstacle that δ n−1 is unknown, we note that correct value of the (n − 2, n − 2) entry of Ĵ n−2 is known; its value can be obtained using (60) but in this case, it can be determined using properties of even and odd functions that its value must be zero. We therefore solve the nonlinear equation where F (δ) is the (n − 2, n − 2) entry of Ĵ n−2 obtained from J n using the above procedure, with δ n−1 = δ.
This equation can be solved using various root-finding methods, such as the secant method. By applying the quadratic formula in solving (61), it can be determined that the solution must lie in (0, 1 / 2]. Choosing initial guesses close to the upper bound of 1 / 2 yields rapid convergence. To improve efficiency, it should be noted that it is not necessary to compute any of the matrices in this algorithm in their entirety to obtain the (n − 2, n − 2) entry of Ĵ n−2 ; only a select few entries from the lower right corner of each matrix are needed. As such, it is possible to solve for δ n−1 in O(1) arithmetic operations, and compute Ĵ n−2 in O(n) operations overall.

Conclusions
We have obtained recurrence relations for generating orthogonal polynomials on the interval (−1, 1) that satisfy the boundary conditions (1) p(1) = 0 and (2) p(−1) = p(1) = 0. These families of orthogonal polynomials can be used to easily implement transformation matrices between physical and frequency space for function spaces of interest for solving PDEs in polar and cylindrical geometries.
While these polynomials are orthogonal with respect to the weight function ω(s) ≡ 1 , it has been shown that they can easily be modified to be orthogonal with respect to the other weight functions. When modified as such to obtain GJPs, recursion coefficients can be obtained with far greater efficiency than by computing the required inner products directly.
Future work includes the development of numerical methods that make use of these families of orthogonal polynomials, or modifications thereof. This includes the adaptation of Krylov subspace spectral methods (Palchak et al. 2015) to polar and cylindrical geometries (Richardson and Lambers 2017).