Sparse spectral and p-finite element methods for partial differential equations on disk slices and trapeziums

Sparse spectral methods for solving partial differential equations have been derived in recent years using hierarchies of classical orthogonal polynomials on intervals, disks, and triangles. In this work we extend this methodology to a hierarchy of non-classical orthogonal polynomials on disk slices (e.g. a half-disk) and trapeziums. This builds on the observation that sparsity is guaranteed due to the boundary being defined by an algebraic curve, and that the entries of partial differential operators can be determined using formulae in terms of (non-classical) univariate orthogonal polynomials. We apply the framework to solving the Poisson, variable coefficient Helmholtz, and Biharmonic equations.


Introduction
This paper develops sparse spectral methods for solving linear partial differential equations on a special class of geometries that includes disk slices and trapeziums. More precisely, we consider the solution of partial differential equations on the domain Ω := {(x, y) ∈ R 2 | α < x < β, γρ(x) < y < δρ(x)} where either of the following conditions hold: Condition 1. ρ is a degree 1 polynomial. Condition 2. ρ is the square root of a non-negative degree ≤ 2 polynomial, −γ = δ > 0.
We show that partial differential equations become sparse linear systems when viewed as acting on expansions involving a family of orthogonal polynomials (OPs) that generalise Jacobi polynomials, mirroring the ultraspherical spectral method for ordinary differential equations [8] and its analogue on the disk [14] and triangle [9,10]. On the disk-slice the family of weights we consider are of the form The corresponding OPs denoted H (a,b,c) n,k (x, y), where n denotes the polynomial degree, and 0 ≤ k ≤ n. We define these to be orthogonalised lexicographically, that is, H (a,b,c) n,k (x, y) = C n,k x n−k y k + (lower order terms) where C n,k = 0 and "lower order terms" includes degree n polynomials of the form x n−j y j where j < k. The precise normalization arises from their definition in terms of onedimensional OPs in Definition 1.
Sparsity comes from expanding the domain and range of an operator using different choices of the parameters a, b and c. Whereas the sparsity pattern and entries derived in [9,10] for equations on the triangle and [14] for equations on the disk results from manipulations of Jacobi polynomials, in the present work we use a more general integration-by-parts argument to deduce the sparsity structure, alongside careful use of the Christoffel-Darboux formula [7, 18.2.2] and quadrature rules to determine the entries. In particular, by exploiting the connection with one-dimensional orthogonal polynomials we can construct discretizations of general partial differential operators of size p(p − 1)/2 × p(p − 1)/2 in O(p 3 ) operations, where p is the total polynomial degree. This compares favourably to O(p 6 ) operations if one proceeds naïvely. Furthermore, we use this framework to derive sparse p-finite element methods that are analogous to those of Beuchler and Schöberl on tetrahedra [1], see also work by Li and Shen [5].
Here is an overview of the paper: Section 2: We present our procedure to gain a (two-parameter) family of 2D orthogonal polynomials (OPs) on the disk-slice domain, by combining 1D OPs on the interval, to form 2D OPs on the disk. Section 3: We demonstrate that these families will lead to sparse operators, including Jacobi operators representing multiplication by x and y, and partial differential operators. We present a method involving use of the Christoffel-Darboux formula [7, 18.2.2] to obtain the recurrence coefficients for the non-classical 1D OPs, allowing us to exactly use the 1D quadrature rules we present to calculate the non-zero entries to the sparse operators. Section 4: We discuss computational issues, in particular, how to realise the results of the preceding sections in practice. We present a method for explicitly deriving the recurrence coefficients of the non-classical 1D OPs we detail in Section 2. We derive a quadrature rule on the disk-slice that can be used to expand a function in the OP basis up to a given order. Further, we implement function evaluation using the coefficients of the expansion of a given function using the Clenshaw algorithm.
Section 5: We demonstrate the proposed technique for solving Poisson, Helmholtz, and Biharmonic equations on the disk-slice.
Appendix A: We use the procedure to construct sparse p-finite element methods. This lays the groundwork for a future hp-finite element method in a disk, where the elements capture the circular geometry precisely.
Appendix B: We discuss extension to the special case of end-disk-slices (e.g., half disks).
Appendix C: We discuss extension to trapezia.

Orthogonal polynomials on the disk-slice and the trapezium
In this section we outline the construction and some basic properties of H (a,b,c,d) n,k (x, y). The symmetry in the weight allows us to express the polynomials in terms of 1D OPs, and deduce certain properties such as recurrence relationships.
be two weight functions on the intervals (α, β) and (γ, δ) respectively, given by: and define the associated inner products by: Denote the three-parameter family of orthonormal polynomials on [α, β] by {R (a,b,c) n }, orthonormal with respect to the inner product defined in (1), and the two-parameter family of orthonormal polynomials on [γ, δ] by {P (a,b) n }, orthonormal with respect to the inner product defined in (2). Definition 2. Define the four-parameter 2D orthogonal polynomials via: , (x, y) ∈ Ω, } are orthogonal with respect to the weight , (x, y) ∈ Ω, assuming that either Condition 1 or Condition 2 with w (a,b) P being an even function (i.e. a = b, and we can hence denote the weight as w (a) where for f, g : Ω → R the inner product is defined as We can see that they are indeed orthogonal using the change of variable t = y ρ(x) , for the following normalisation: For the disk-slice, the weight Note here we can simply remove the need for including a fourth parameter d. The 2D OPs orthogonal with respect to the weight above on the disk-slice Ω are then given by: , (x, y) ∈ Ω In this case the weight w P (x) is an ultraspherical weight, and the corresponding OPs are the normalized Jacobi polynomials {P (b,b) n }, while the weight w R (x) is non-classical (it is in fact semi-classical, and is equivalent to a generalized Jacobi weight [6, §5]).

Jacobi matrices
We can express the three-term recurrences associated with R (a,b,c) n and P Of course, for the disk-slice case, we have that c = d and γ (c,c) n = 0 ∀n = 0, 1, 2, . . . . We can use (7) and (8) to determine the 2D recurrences for H (a,b,c,d) n,k (x, y). Importantly, we can deduce sparsity in the recurrence relationships: n,k (x, y) satisfy the following 3-term recurrences: Proof. The 3-term recurrence for multiplication by x follows from equation (7). For the recurrence for multiplication by y, since {H Recall from equation (5) . Then for m = 0, . . . , n + 1, j = 0, . . . , m, using the change of variable t = y ρ(x) : where, by orthogonality, Finally, if Condition 1 holds we have that If Condition 2 holds we have that γ Three-term recurrences lead to Jacobi operators that correspond to multiplication by x and y. Define, for n = 0, 1, 2, . . . : as the Jacobi matrices corresponding to The matrices J are banded-block-banded matrices: and sub-block-bandwidths (λ, µ) if all blocks A i,j are banded with bandwidths (λ, µ). A matrix where the block-bandwidths and sub-block-bandwidths are small compared to the dimensions is referred to as a banded-block-banded matrix.
are block-tridiagonal (block-bandwidths (1, 1)): where the blocks themselves are diagonal for J (a,b,c,d) x (sub-block-bandwidths (0, 0)), Note that the sparsity of the Jacobi matrices (in particular the sparsity of the sub-blocks) comes from the natural sparsity of the three-term recurrences of the 1D OPs, meaning that the sparsity is not limited to the specific disk-slice case.

Building the OPs
We can combine each system in (9) into a block-tridiagonal system: , and for each n = 0, 1, 2 . . . , For each n = 0, 1, 2 . . . let D n be any matrix that is a left inverse of A n , i.e. such that D n A n = I n+2 . Multiplying our system by the preconditioner matrix that is given by the block diagonal matrix of the D n 's, we obtain a lower triangular system [4, p78], which can be expanded to obtain the recurrence: Note that we can define an explicit D n as follows: c basis conversion where It follows that we can apply D n in O(n) complexity, and thereby calculate H n,k,8 ≡ 0 for any n, k.

Sparse partial differential operators
In this section, we concentrate on the disk-slice case, and simply note that similar arguments apply for the trapezium case. Recall that, for the disk-slice, The 2D OPs on the disk-slice Ω, orthogonal with respect to the weight are then given by: Denote the weighted OPs by and recall that a function f (x, y) defined on Ω is approximated by its expansion f (x, according to: The incrementing and decrementing of parameters as seen here is analogous to other well known orthogonal polynomial families' derivatives, for example the Jacobi polynomials on the interval, as seen in the DLMF [7, (18.9.3)], and on the triangle [9].
has block-bandwidths (3, −1), and sub-block-bandwidths (2, 0). Proof. First, note that: We proceed with the case for the operator D . The coefficients of the expansion are then the entries of the relevant operator matrix. We can use an integrationby-parts argument to show that the only non-zero coefficient of this expansion is when Then, using the change of variable t = y ρ(x) , we have that P is a polynomial of degree 2c and vanishes at the limits of the integral for positive parameter c, we have that which is zero for j < k − 1 by orthogonality. Further, when j = k − 1, we have that showing that the only possible non-zero coefficient is when We next consider the case for the operator D . The coefficients of the expansion are then the entries of the relevant operator matrix. As before, we can use an integration-by-parts argument to show that the only non-zero coefficients of this expansion are when m = n − 1, n − 2, n − 3, j = k, k − 1, k − 2 and 0 ≤ j ≤ m. First, note that We will first show that the second factor of each term in (13) are zero for j < k − 2 and also for j = k − 1. To this end, observe that, for any integer c, P (c,c) (−t) = (−1) k P (c,c) (t) and so P (c,c) k is an even polynomial for even k, and an odd polynomial for odd k. Thus, P is an odd polynomial for any k. Hence is zero for j < k − 2 by orthogonality, and is zero for j = k − 1 due to symmetry over the domain. Moreover, t P (c,c) k (t) P (c+1,c+1) j (t) is also an odd polynomial for any k and so is zero for j = k − 1 due to symmetry over the domain, and which is zero for j < k − 2 by orthogonality. Thus, (13) is zero for j / ∈ {k − 2, k}.
Now, using (10), integration-by-parts, and noting that the weight w (a,b,2c) R is a polynomial degree a + b + 2c and vanishes at the limits of the integral for positive parameters a, b, c, we have that By recalling (12) and noting that j − k is even by the earlier argument, we can see ρ ρ w Hence, by orthogonality, each term in (14) is is zero for m − j + 3 + j − k < n − k ⇐⇒ m < n − 3. Finally, which is also zero for m < n − 3. Thus showing that the only possible non-zero coefficients c x m,j are when m = n − 3, . . . , n and j = k − 2, k (j ≤ m).
We can gain the non-zero entries of the weighted differential operators similarly, by noting that for the disk-slice and also that There exist conversion matrix operators that increment/decrement the parameters, transforming the OPs from one (weighted or non-weighted) parameter space to another. for conversion between weighted spaces, according to: Lemma 2. The operator matrices in Definition 5 are sparse, with banded-block-banded structure. More specifically: • T (a,b,c)→(a+1,b+1,c) has block-bandwidth (0, 2), with diagonal blocks.
• T = ω Since w (c) P is a polynomial degree 2c, we have that the above is then zero for j < k − 2c. Further, since w (ã,b, 2c+j−k) R is a polynomial of degreeã +b + 2c + j − k, we have that the above is zero for m − j +ã +b + 2c + j − k < n − k ⇐⇒ m < n −ã −b − 2c.
The sparsity argument for the weighted parameter transformation operators follows similarly. General linear partial differential operators with polynomial variable coefficients can be constructed by composing the sparse representations for partial derivatives, conversion between bases, and Jacobi operators. As a canonical example, we can obtain the matrix operator for the Laplacian ∆, that will take us from coefficients for expansion in the weighted space to coefficients in the non-weighted space H (1,1,1) (x, y). Note that this construction will ensure the imposition of the Dirichlet zero boundary conditions on Ω. The matrix operator for the Laplacian we denote ∆ Importantly, this operator will have banded-block-banded structure, and hence will be sparse, as seen in Figure 2.
Another important example is the Biharmonic operator ∆ 2 , where we assume zero Dirichlet and Neumann conditions. To construct this operator, we first note that we can obtain the matrix operator for the Laplacian ∆ that will take us from coefficients for expansion in the space H (0,0,0) (x, y) to coefficients in the space H (2,2) (x, y). We denote this matrix operator that acts on the coefficients vector as ∆ (0,0,0)→ (2,2,2) , and is given by Further, we can represent the Laplacian as a map from coefficients in the space W (2,2) to coefficients in the space H (0,0,0) . Note that a function expanded in the W (2,2) basis will satisfy both zero Dirichlet and Neumann boundary conditions on Ω. We denote this matrix operator as ∆ (2,2,2)→(0,0,0) W , and is given by We can then construct a matrix operator for ∆ 2 that will take coefficients in the space W (2,2,2) to coefficients in the space H (2,2,2) . Note that any function expanded in the W (2,2,2) basis will satisfy both zero Dirichlet and zero Neumann boundary conditions on Ω. The matrix operator for the Biharmonic operator is then given by

Computational aspects
In this section we discuss how to take advantage of the proposed basis and sparsity structure in partial differential operators in practical computational applications.

Constructing
It is possible to obtain the recurrence coefficients for the {R (a,b,c) n } OPs in (7), by careful application of the Christoffel-Darboux formula [7, 18.2.12]. We explain the process here for the disk-slice case, however we note that a similar but simpler argument holds for the trapezium case. We thus first need to define a new set of 'interim' 1D OPs.
be a weight function on the interval (α, β), and define the associated inner product by: Proof. Fix n, m ∈ {0, 1, . . . } and without loss of generality, assume m ≤ n. First recall that and define Now, by the Christoffel-Darboux formula [7, 18.2.12], we have that for any x, y ∈ R, Then, using (22) and (24), showing that the RHS and LHS of (20) are equivalent. Further, } are given by: (1)β Let f (x, y) be a polynomial on Ω. The quadrature rule is then where M = 1 2 (N + 1) 2 , and the quadrature rule is exact if f (x, y) is a polynomial of degree ≤ N .
Proof. We will use the substitution that First, note that, for (x, y) ∈ Ω, so that y → f e (x, y) for fixed x is an even function, and y → f o (x, y) for fixed x is an odd function. Note that if f is a polynomial, then f e (s, ρ(s)t) is a polynomial in s ∈ [α, β] for fixed t.

Now, we have that
Suppose f is a polynomial in x and y of degree N , and hence that f e is a degree ≤ N polynomial. First, note that the degree of the polynomial given by x → f e (x, y) for fixed y is ≤ N and the degree of the polynomial given by y → f e (x, y) for fixed x is ≤ N . Also note that s → f e s, ρ(s)t for fixed t is then a degree N polynomial (since ρ is a degree 1 polynomial). Hence, we achieve equality at ( ) if 2M 2 − 1 ≥ N and we achieve equality at ( ) if also 2M 1 − 1 ≥ N .
Next, note that  Hence, for a polynomial f in x and y of degree N , where M = 1 2 (N + 1) 2 .

Examples on the disk-slice with zero Dirichlet conditions
We now demonstrate how the sparse linear systems constructed as above can be used to efficiently solve PDEs with zero Dirichlet conditions. We consider Poisson, inhomogeneous variable coefficient Helmholtz equation and the Biharmonic equation, demonstrating the versatility of the approach.

Poisson
The Poisson equation is the classic problem of finding u(x, y) given a function f (x, y) such that: noting the imposition of zero Dirichlet boundary conditions on u.
We can tackle the problem as follows. Denote the coefficient vector for expansion of u in the W (1,1,1) OP basis up to degree N by u, and the coefficient vector for expansion of f in the H (1,1,1) OP basis up to degree N by f . Since f is known, we can obtain f using the quadrature rule above. In matrix-vector notation, our system hence becomes: which can be solved to find u. In Figure 3 we see the solution to the Poisson equation with zero boundary conditions given in (37) in the disk-slice Ω. In Figure 3 we also show the norms of each block of calculated coefficients of the approximation for four right-hand sides of the Poisson equation with N = 200, that is, 20,301 unknowns. The rate of decay in the coefficients is a proxy for the rate of convergence of the computed solution: as typical of spectral methods, we expect the numerical scheme to converge at the same rate as the coefficients decay. We see that we achieve algebraic convergence for the first three examples, noting that for right hand-sides that vanish at the corners of our disk-slice (x ∈ {α, β}, y = ±ρ(x)) we observe faster convergence.
In Figure 4 we see an example where the solution calculated to the Poisson equation is shown together with a plot of the exact solution and the error. The example was chosen so that the exact solution was u(x, y) = W (1,1,1) (x, y)y 3 exp(x), and thus the RHS function f would be f (x, y) = ∆ W (1,1,1) (x, y)y 3 exp(x) . We see that the computed solution is almost exact.
where k ∈ R, noting the imposition of zero Dirichlet boundary conditions on u.
We can tackle the problem as follows. Denote the coefficient vector for expansion of u in the W (1,1,1) OP basis up to degree N by u, and the coefficient vector for expansion of f in the H (1,1,1) OP basis up to degree N by f . Since f is known, we can obtain the coefficients f using the quadrature rule above. We can obtain the matrix operator for the variable-coefficient function v(x, y) by using the Clenshaw algorithm with matrix inputs as the Jacobi matrices J which can be solved to find u. We can see the sparsity and structure of this matrix system in Figure 2 with v(x, y) = xy 2 as an example. In Figure 5 we see the solution to the inhomogeneous variable-coefficient Helmholtz equation with zero boundary conditions given in (38) in the half-disk Ω, with k = 100, v(x, y) = 1 − (3(x − 1) 2 + 5y 2 ) and f (x, y) = x(1 − x 2 − y 2 )e x . In Figure 5 we also show the norms of each block of calculated coefficients of the approximation for four right-hand sides of the inhomogeneous variablecoefficient Helmholtz equation with k = 20 and v(x, y) = 1 − (3(x − 1) 2 + 5y 2 ) using N = 200, that is, 20,301 unknowns. The rate of decay in the coefficients is a proxy for the rate of convergence of the computed solution. We see that we achieve algebraic convergence for the first three examples, noting that for right hand sides that vanish at the corners of our disk-slice (x ∈ {α, β}, y = ±ρ(x)) we see faster convergence.

Biharmonic equation
Find u(x, y) given a function f (x, y) such that: where ∆ 2 is the Biharmonic operator, noting the imposition of zero Dirichlet and Neumann boundary conditions on u. In Figure 6 we see the solution to the Biharmonic equation (39) in the disk-slice Ω. In Figure 6 we also show the norms of each block of calculated coefficients of the approximation for four right-hand sides of the biharmonic equation with N = 200, that is, 20,301 unknowns. We see that we achieve algebraic convergence for the first three examples, noting that for right hand sides that vanish at the corners of our disk-slice (x ∈ {α, β}, y = ±ρ(x)) we see faster convergence.

Conclusions
We have shown that bivariate orthogonal polynomials can lead to sparse discretizations of general linear PDEs on specific domains whose boundary is specified by an algebraic curve-notably here the disk-slice-with Dirichlet boundary conditions. This work extends the triangle case [1,5,10] to non-classical geometries, and forms a building block in developing an hp−finite element method to solve PDEs on other polygonal domains by using suitably shaped elements, for example, by dividing the disk into disk slice elements. This work serves as a stepping stone to constructing similar methods to solve partial differential equations on 3D sub-domains of the sphere, such as spherical caps and spherical triangles. In particular, orthogonal polynomials (OPs) in cartesian coordinates (x, y, and z) on a half-sphere can be represented using two families of OPs on the half-disk, see [11, Theorem 3.1] for a similar construction of OPs on an arc in 2D, and it is clear from the construction in this paper that discretizations of spherical gradients and Laplacian's are sparse on half-spheres and other suitable sub-components of the sphere. The resulting sparsity in high-polynomial degree discretizations presents an attractive alternative to methods based on bijective mappings (e.g., [2,12,3]). Constructing these sparse spectral methods for surface PDEs on half-spheres, spherical caps, and spherical triangles is future work, and has applications in weather prediction [13]. Other extensions include a full hp-finite element method on sections of a disk, which has applications in turbulent pipe flow.

A P-finite element methods using sparse operators
We follow the method of [1] to construct a sparse p-finite element method in terms of the operators constructed above, with the benefit of ensuring that the resulting discretisation is symmetric. Consider the 2D Dirichlet problem on a domain Ω: In general, we would let T be the set of elements τ that make up our finite element discretisation of the domain, where each τ is a trapezium or disk slice for example.
In this section, we limit our discretisation to a single element, that is we let τ = Ω for a diskslice domain. We can choose our finite dimensional space V p = {v p ∈ V | deg (v p | τ ) ≤ p} for some p ∈ N.
We seek u p ∈ V p s.t.
Define Λ (a,b,c) : = H (a,b,c) , H (a,b,c) W (a,b,c) where W (a,b,c) is the weight with which the OPs in H (a,b,c) are orthogonal with respect to. Note that due to orthogonality this is a diagonal matrix. We can choose a basis for V p by using the weighted orthogonal polynomials on τ with parameters a = b = 1: 0 (x, y) W (1,1,1) 1 (x, y) W (1,1,1) 2 (x, y) . . . Note that since we have sparse operator matrices for partial derivatives and basis-transform, we obtain a symmetric sparse (element) stiffness matrix, as well as a sparse operator matrix for calculating the load vector (rhs).

B End-Disk-Slice
The work in this paper on the disk-slice can be easily transferred to the special-case domain of the end-disk-slice , such as half disks, by which we mean Ω := {(x, y) ∈ R 2 | α < x < β, γρ(x) < y < δρ(x)}