A sparse spectral method on triangles

TL;DR

This paper develops a spectral method using bivariate orthogonal polynomials on triangles, enabling efficient, sparse discretizations for solving large-scale linear PDEs with high polynomial degrees.

Contribution

It introduces a spectral method on triangles leveraging analogues of univariate tools, allowing for sparse, high-degree polynomial discretizations of PDEs.

Findings

Enables solving PDEs with polynomial degrees in the thousands.

Achieves sparse discretizations with millions of degrees of freedom.

Provides practical algorithms for spectral methods on triangular domains.

Abstract

In this paper, we demonstrate that many of the computational tools for univariate orthogonal polynomials have analogues for a family of bivariate orthogonal polynomials on the triangle, including Clenshaw's algorithm and sparse differentiation operators. This allows us to derive a practical spectral method for solving linear partial differential equations on triangles with sparse discretizations. We can thereby rapidly solve partial differential equations using polynomials with degrees in the thousands, resulting in sparse discretizations with as many as several million degrees of freedom.

Equations

$${\Delta }^{2}u$$

$$u_y$$

$$\Delta u+V(x,y)u$$

$$P_{n,k}(x,y)$$

$${\iint }_T{P}_{n,k}(x,y)P_{m,\ell }(x,y)\mathrm{d}x\mathrm{d}y=\{\begin{array}{ll}\frac{1}{{\pi }_{n,k}},& (n,k)=(m,\ell ),\\ 0,& (n,k)\ne (m,\ell ),\end{array}$$

$$f(x,y)=\sum _{n=0}^{\infty }\sum _{k=0}^n{f}_{n,k}{P}_{n,k}(x,y),f_{n,k}={\pi }_{n,k}{\iint }_{T}f(x,y)P_{n,k}(x,y)\mathrm{d}x\mathrm{d}y,$$

$$f(x,y)\approx \sum _{n=0}^N\sum _{k=0}^n{a}_{n,k}{P}_{n,k}(x,y)$$

$$P_{n,k}^{(a,b,c)}(x,y)={\tilde{P}}_{n-k}^{(2k+b+c+1,a)}(x)(1-x{)}^k{\tilde{P}}_k^{(c,b)}\left(\frac{y}{1-x}\right),n\ge k\ge 0,$$

$$u=\sum _{n=0}^N\sum _{k=0}^n{a}_{n,k}{P}_{n,k}\Rightarrow \frac{\partial u}{\partial y}=\sum _{n=0}^{N-1}\sum _{k=0}^n(k+1)a_{n+1,k+1}{P}_{n,k}^{(0,1,1)}.$$

$$p_{0,0}(x,y),p_{1,0}(x,y),p_{1,1}(x,y),p_{2,0}(x,y),p_{2,1}(x,y),p_{2,2}(x,y),\dots ,$$

$${\iint }_{\Omega }w(x,y)p_{n,k}(x,y)p_{m,\ell }(x,y)\mathrm{d}x\mathrm{d}y=\{\begin{array}{ll}{d}_{n,k},& (n,k)=(m,\ell ),\\ 0,& (n,k)\ne (m,\ell ),\end{array}$$

$${\mathbb{P}}_n(x,y)=\left(\begin{array}{c}{p}_{n,0}(x,y)\\ ⋮\\ p_{n,n}(x,y)\end{array}\right).$$

$${\iint }_{\Omega }w(x,y){\mathbb{P}}_m(x,y){\mathbb{P}}_n(x,y{)}^{\top }\mathrm{d}x\mathrm{d}y=\{\begin{array}{ll}{D}_n,& m=n,\\ 0,& m\ne n,\end{array}$$

$$\mathbf{P}(x,y)={\left({\mathbb{P}}_0(x,y\right),{\mathbb{P}}_1(x,y),\dots )}^{\top }$$

$$f(x,y)$$

$$f_{n,k}=\frac{1}{d_{n,k}}{\iint }_{\Omega }w(x,y)f(x,y)p_{n,k}(x,y)\mathrm{d}x\mathrm{d}y,n\ge k\ge 0,$$

$$b_k{p}_{k+1}(x)+a_k{p}_k(x)+c_{k-1}{p}_{k-1}(x)=x{p}_k(x)$$

$$J\left(\begin{array}{c}{p}_0(x)\\ p_1(x)\\ ⋮\end{array}\right)=x\left(\begin{array}{c}{p}_0(x)\\ p_1(x)\\ ⋮\end{array}\right),J=\left(\begin{array}{cccc}{a}_0& b_0& & \\ c_0& a_1& b_1& \\ & c_1& a_2& \ddots \\ & & \ddots & \ddots \end{array}\right).$$

$$xf(x)=x{\mathbf{P}}(x,y)^{\top}\mbox{\boldmath$f$\unboldmath}={\mathbf{P}}(x,y)^{\top}J^{\top}\mbox{\boldmath$f$\unboldmath}.$$

$$J_x\mathbf{P}(x,y)=x\mathbf{P}(x,y),J_y\mathbf{P}(x,y)=y\mathbf{P}(x,y).$$

$$J_x=\left(\begin{array}{cccc}{A}_0^x& B_0^x& & \\ C_0^x& A_1^x& B_1^x& \\ & C_1^x& A_2^x& \ddots \\ & & \ddots & \ddots \end{array}\right),J_y=\left(\begin{array}{cccc}{A}_0^y& B_0^y& & \\ C_0^y& A_1^y& B_1^y& \\ & C_1^y& A_2^y& \ddots \\ & & \ddots & \ddots \end{array}\right),$$

$$x{\mathbf{P}}(x,y)^{\top}\mbox{\boldmath$f$\unboldmath}={\mathbf{P}}(x,y)^{\top}J_{x}^{\top}\mbox{\boldmath$f$\unboldmath}\quad\hbox{and}\quad y{\mathbf{P}}(x,y)^{\top}\mbox{\boldmath$f$\unboldmath}={\mathbf{P}}(x,y)^{\top}J_{y}^{\top}\mbox{\boldmath$f$\unboldmath}.$$

$$L_N(x_{\ast })\left(\begin{array}{c}{p}_0(x_{\ast })\\ p_1(x_{\ast })\\ p_2(x_{\ast })\\ ⋮\\ p_N(x_{\ast })\end{array}\right)=\left(\begin{array}{ccccc}1& & & & \\ a_0-x_{\ast }& b_0& & & \\ c_0& a_1-x_{\ast }& b_1& & \\ & \ddots & \ddots & \ddots & \\ & & c_{N-2}& a_{N-1}-x_{\ast }& b_{N-1}\end{array}\right)\left(\begin{array}{c}{p}_0(x_{\ast })\\ p_1(x_{\ast })\\ p_2(x_{\ast })\\ ⋮\\ p_N(x_{\ast })\end{array}\right)=e_0,$$

$$f(x_{\ast })={\left(p_0(x_{\ast }\right),\dots ,p_N(x_{\ast }))}^{\top }\mathbf{a}=e_0^{\top }\left({\left(L_N(x_{\ast })\right)}^{-\top }\mathbf{a}\right),\mathbf{a}=\left(\begin{array}{c}{a}_0\\ ⋮\\ a_N\end{array}\right).$$

$$L_N(x_{\ast },y_{\ast })\mathbf{P}(x_{\ast },y_{\ast })=\left(\begin{array}{cccc}1& & & \\ A_0^x-x_{\ast }{I}_1& B_0^x& & \\ A_0^y-y_{\ast }{I}_1& B_0^y& & \\ C_0^x& A_1^x-x_{\ast }{I}_2& B_1^x& \\ C_0^y& A_1^y-y_{\ast }{I}_2& B_1^y& \\ & \ddots & \ddots & \ddots \end{array}\right)\mathbf{P}(x_{\ast },y_{\ast })=\left(\begin{array}{c}1\\ 0_{1\times 1}\\ 0_{1\times 1}\\ 0_{2\times 1}\\ 0_{2\times 1}\\ ⋮\end{array}\right),$$

$$B_n=\left(\begin{array}{c}{B}_n^x\\ B_n^y\end{array}\right)\in {\mathbb{R}}^{(2n+2)\times (n+2)}$$

$${\tilde{L}}_N(x_{\ast },y_{\ast })\mathbf{P}(x_{\ast },y_{\ast })=\left(\begin{array}{c}1\\ 0_{1\times 1}\\ 0_{2\times 1}\\ ⋮\end{array}\right),{\tilde{L}}_N(x_{\ast },y_{\ast })=\left(\begin{array}{cccc}1& & & \\ & B_0^{+}& & \\ & & B_1^{+}& \\ & & & \ddots \end{array}\right)L_N(x_{\ast },y_{\ast }).$$

$$f(x_{*},y_{*})=\mathbf{P}(x_{*},y_{*})^{\top}\mbox{\boldmath$a$\unboldmath}={\mbox{\boldmath$e$\unboldmath}}_{0}^{\top}\!\left(\left(\tilde{L}_{N}(x_{*},y_{*})\right)^{-\top}\!\mbox{\boldmath$a$\unboldmath}\right).$$

$$\mbox{\boldmath$g$\unboldmath}=M_{q}\mbox{\boldmath$f$\unboldmath},\qquad M_{q}=q(J_{x}^{\top},J_{y}^{\top}),$$

$$q(J_x^{\top },J_y^{\top })=\sum _{n=0}^d\sum _{k=0}^n{c}_{nk}(J_x^{\top }{)}^{n-k}(J_y^{\top }{)}^k,q(x,y)=\sum _{n=0}^d\sum _{k=0}^n{c}_{nk}{x}^{n-k}{y}^k.$$

$$M_q=(e_0\otimes \mathcal{I})(L^{-\top }\mathbf{q}),L=\left(\begin{array}{cccc}{I}_1\otimes \mathcal{I}& & & \\ A_0^x\otimes \mathcal{I}-I_1\otimes J_x& B_0^x\otimes \mathcal{I}& & \\ A_0^y\otimes \mathcal{I}-I_1\otimes J_y& B_0^y\otimes \mathcal{I}& & \\ C_0^x\otimes \mathcal{I}& A_1^x\otimes \mathcal{I}-I_2\otimes J_x& B_1^x\otimes \mathcal{I}& \\ C_0^y\otimes \mathcal{I}& A_1^y\otimes \mathcal{I}-I_2\otimes J_y& B_1^y\otimes \mathcal{I}& \\ & \ddots & \ddots & \ddots \end{array}\right),$$

Full text

{AMS}

33D50, 65N35

1 Introduction

Univariate orthogonal polynomials are fundamental in applied and computational mathematics. They are used for the development of quadrature rules [8], spectral theory of Jacobi operators [24], eigenvalue statistics of random matrices [5], computational approximation theory [26], and to derive spectral methods for the numerical solution of differential equations [2, 3, 15, 18, 27, 28]. On the contrary, multivariate orthogonal polynomials currently have a more limited impact in applications and computational methods, though it is an active research area with a promising future.

To demonstrate the potential practical importance of multivariate orthogonal polynomials, we show that many computational tools for univariate orthogonal polynomials can be generalized to a family of bivariate orthogonal polynomials on a triangle. These tools allow us to derive a sparse spectral method for solving general linear partial differential equations (PDEs) with Dirichlet and Neumann conditions on triangles. While the techniques are general, we demonstrate the method on the following PDEs:

[TABLE]

Since triangles can be mapped to each other by affine translations, and polynomials remain polynomial, we can consider a single reference triangle, without loss of generality. Throughout this paper, we select the reference triangle to be the unit simplex: a right-angled triangle of unit height and width, i.e., $T=\{(x,y):0<x<1,0<y<1-x\}$.

There are several different families of bivariate orthogonal polynomials on $T$ [7]. Here, we consider a family that is built from univariate orthogonal polynomials [11]:

[TABLE]

where $\tilde{P}_{k}^{(a,b)}(x)$ denotes the degree $k$ shifted Jacobi polynomial on $[0,1]$ with parameters $(a,b)$.111In particular, $\tilde{P}_{k}^{(a,b)}(x)=P_{k}^{(a,b)}(2x-1)$, where $P_{k}^{(a,b)}$ is the degree $k$ Jacobi polynomial on $[-1,1]$ with parameters $(a,b)$. The polynomials in Eq. 1 are one possible generalization on triangles of the Legendre polynomials [13, Tab. 18.3.1]. In particular, the polynomials satisfy three-term recurrence relations (see Eq. 7) and are orthogonal with respect to the standard $L^{2}$ inner-product on $T$:

[TABLE]

where $\pi_{n,k}=2(2k+1)(n+1)$. They provide a well-conditioned basis to represent integrable functions $f\in L^{2}(T)$ as a series expansion,

[TABLE]

where the first equality above should be understood in the $L^{2}$-sense. In order to do efficient computations with functions defined on a triangle, it is important to be able to rapidly compute expansion coefficients of $f(x,y)$ so that

[TABLE]

for a selected integer $N$. Recently, Slevinsky developed and implemented a fast backward stable algorithm for precisely this task [19, 20], accompanied with an optimized multithreaded open-source C library [21], allowing expansions to be computationally feasible for relatively large $N$. This has greatly improved the practicality of spectral methods for triangular domains.

The use of bivariate orthogonal polynomials on triangles has a long history in the spectral element method and $p$-finite element method ($p$-FEM) literature [10], going back to Dubiner [6]. The polynomials in Eq. 1 lead to highly structured $p$-FEM discretization matrices for PDEs of the form $\mathcal{L}u=-\nabla\cdot(A(x,y)\nabla u)$, and when $A(x,y)$ is a constant one can derive sparse discretizations that can be generated in optimal complexity [1]. The present work can be viewed as a generalization of this construction to strong formulations of PDEs that are not necessarily elliptic. Moreover, the properties of bivariate orthogonal polynomials allows us to retain sparsity for high differential order and variable coefficient PDEs (see Section 4.1.2).

Our main idea is to exploit a hierarchy of sparse recurrence relations [16] that hold between the polynomials in Eq. 1 and the so-called Jacobi polynomials on the triangle [7, 11]:222The polynomials $P_{n,k}^{(a,b,c)}$ for $a,b,c>-1$ satisfy $\iint_{T}P_{n,k}^{(a,b,c)}(x,y)P_{m,\ell}^{(a,b,c)}(x,y)x^{a}y^{b}(1-x-y)^{c}\text{d}x\text{d}y=0$ if $n\neq m$ or $k\neq\ell$.

[TABLE]

where $a,b,c>-1$. In a manner that is analogous to the ultraspherical spectral method [15, 25], we represent the action of partial derivatives by representing the domain and range as vectors of coefficients in different bases so that the matrix representation is sparse. For example, while $\smash{\tfrac{\partial}{\partial y}P_{n,k}}$ for $k\geq 1$ cannot be written as a sparse vector of $P_{n,k}$ coefficients, we have $\smash{\frac{\partial}{\partial y}P_{n,k}=(k+1)P_{n-1,k-1}^{(0,1,1)}}$, (see Corollary A.1). This means that the first partial derivative with respect to $y$ has a sparse matrix representation if the range is represented as a vector of $P_{n,k}^{(0,1,1)}$ coefficients. This can be summarized as

[TABLE]

Moreover, these sparse recurrence relationships form a hierarchy, in the sense that $\smash{\tfrac{\partial^{s}}{\partial y^{s}}}$ has a sparse representation if the range is represented as a vector of $\smash{P_{n,k}^{(0,s,s)}}$ coefficients, for any $s\geq 0$. Similar, but slightly more complicated, sparse recurrence relations hold for $\smash{\tfrac{\partial^{s}}{\partial x^{s}}P_{n,k}}$ when the range is represented as vectors in $\smash{P_{n,k}^{(s,0,s)}}$ coefficients for any $s\geq 0$ (see Section 3.4).

One is also able to combine sparse representations to discretize linear PDEs. For example, the Laplacian operator $\Delta u=u_{xx}+u_{yy}$ can be represented by a sparse matrix if the range is selected to be a vector of $\smash{P_{n,k}^{(2,2,2)}}$ coefficients while the domain is a vector of $P_{n,k}$ coefficients. This is because there exist sparse conversion relationships for converting between certain $\smash{P_{n,k}^{(a,b,c)}}$ bases (see Section 3.1). Figure 1 illustrates a typical schema that illustrates how sparse recurrences are combined. In the language of finite-element methods, the test and trial spaces are different with a sparse embedding of the trial space in the test space.

The paper is organized as follows. In Section 2, we establish some general computational tools for bivariate orthogonal polynomials such as Jacobi operators and the bivariate Clenshaw algorithm. In Section 3, we specialize to Eq. 2, where the additional structure allows us to achieve a more efficient Clenshaw algorithm. In Section 4, we employ weighted Jacobi polynomials on the triangle to solve PDEs such as a variable coefficient Helmholtz equation and a biharmonic equation with zero Dirichlet conditions. In Section 5, we extend the ideas to solve linear PDEs with nonzero Dirichlet conditions, and in Section 6 we demonstrate that the framework easily generalizes to systems of PDEs so that it can be used to solve the Helmholtz equation in a polygonal domain.

The appendices contain relationships and additional formulae about orthogonal polynomials on the triangle. Our spectral method depends on explicit rational recurrence relationships that the polynomials $P_{n,k}^{(a,b,c)}(x,y)$ satisfy for differentiation, weighted differentiation, and conversion, which we detail in Appendix A. Tackling Dirichlet conditions requires a modification of the basis to enable sparse restriction operators, which we define as $Q_{n,k}^{(a,b,c)}$ in Appendix B. These also have explicit rational recurrence relationships for differentiation and conversion, which we derive in Appendix C.

2 Computations with bivariate orthogonal polynomials

In this section, we derive several computational tools for bivariate orthogonal polynomials such as the Jacobi operators, Clenshaw’s algorithm, and multiplication operators. Later, in Section 3, we specialize these tools to the Jacobi polynomials on the triangle (see Eq. 2).

Consider a sequence of bivariate polynomials

[TABLE]

where $\left\{p_{n,k}\right\}_{0\leq k\leq n\leq N}$ is a basis for the space of bivariate polynomials of total degree $\leq N$,333We say that a bivariate polynomial $q(x,y)$ is of total degree $\leq N$ if $q(x,y)=\sum_{n=0}^{N}\sum_{k=0}^{n}b_{n,k}x^{k}y^{n-k}$ for some coefficients $b_{n,k}$. for any integer $N$. We say that such a sequence is orthogonal with respect to a nonnegative weight function $w(x,y)$ on $\Omega\subset\mathbb{R}^{2}$ if

[TABLE]

where $d_{n,k}$ are positive numbers.

It is notationally convenient to write the bivariate polynomials of the same total degree as a single vector-valued polynomial [7] as follows:

[TABLE]

One can then state the orthogonality condition in Eq. 3 more succinctly as

[TABLE]

where $D_{n}$ is the $(n+1)\times(n+1)$ diagonal matrix with entries $d_{n,k}$ for $0\leq k\leq n$, $\mathbf{0}$ is a matrix of all zeros of the appropriate size, and ${\mathbb{P}}_{n}(x,y)^{\top}$ denotes the transpose of ${\mathbb{P}}_{n}(x,y)$. The sequence of bivariate polynomials are normalized (orthonormal) if $D_{n}$ is the identity matrix for all $n\geq 0$. We also use the notation

[TABLE]

to encode all of the polynomials as a single infinite vector.

2.1 Bivariate function approximation

A sequence of bivariate orthogonal polynomials on $\Omega\subset\mathbb{R}^{2}$ can be used to approximate functions that are square integrable with respect to the associated weight function $w(x,y)$ on $\Omega$. For example, provided $\iint_{\Omega}w(x,y)\left|f(x,y)\right|^{2}{\rm d}x{\rm d}y<\infty$, we can write

[TABLE]

where $\boldsymbol{f}_{n}\!=\!{{\left({f_{n,0},\ldots,f_{n,n}}\right)}^{\top}}$​ and \mbox{\boldmathf\unboldmath}\!=\!{{\left({\boldsymbol{f}_{0},\boldsymbol{f}_{1},\ldots}\right)}^{\top}}​ are the coefficients of the expansion. Here, the first equality in Eq. 5 is understood in the sense that the difference between the left- and right-hand side is zero in the norm associated to the inner-product.

The expansion coefficients in Eq. 5 are defined by the following integrals:

[TABLE]

where $d_{n,k}$ is the orthogonality constant in Eq. 3. In practice, it is usually desirable for the expansion coefficients to rapidly decay, i.e., $\|\boldsymbol{f}_{0}\|,\|\boldsymbol{f}_{1}\|,\ldots$ is a rapidly decaying sequence.

2.2 Jacobi operators

In the theory of univariate orthogonal polynomial an important object is the Jacobi operator, which is a self-adjoint linear operator given by a tridiagonal matrix [24]. It is closely related to the fact that a sequence of univariate orthogonal polynomials satisfies a three-term recurrence. For example, if $p_{0},p_{1},\ldots,$ is a sequence of univariate orthogonal polynomials, then

[TABLE]

for $k\geq 1$ [23, Thm. 3.2.1] and

[TABLE]

The Jacobi operator associated with $p_{0},p_{1},\ldots,$ is the symmetric tridiagonal matrix obtained by a diagonal similarity transform of $J$ [24]. This diagonal similarity transform corresponds precisely to the normalization factors required to orthonormalize the sequence of univariate orthogonal polynomials. The transformation is possible provided $0<b_{k}^{-1}c_{k}<\infty$ for all $k$. In particular, if $\{p_{k}(x)\}_{k\geq 0}$ are orthonormal, then $J$ is a symmetric tridiagonal matrix.

A related fact that is important for designing spectral methods is that $J^{\top}$ can be interpreted as the “multiplication-by-$x$” operator. That is, if f(x)={\mathbf{P}}(x,y)^{\top}\mbox{\boldmathf\unboldmath} we have

[TABLE]

In other words, J^{\top}\mbox{\boldmathf\unboldmath} gives the coefficients of $xf(x)$.

The analogue for bivariate orthogonal polynomials is a pair of commuting operators $J_{x}$ and $J_{y}$ [7, §3.4], which satisfy

[TABLE]

Here, $J_{x}$ and $J_{y}$ are block tridiagonal operators so that

[TABLE]

where $A_{n}^{x},A_{n}^{y}\in{\mathbb{R}}^{(n+1)\times(n+1)}$ , $B_{n}^{x},B_{n}^{y}\in{\mathbb{R}}^{(n+1)\times(n+2)}$, and $C_{n}^{x},C_{n}^{y}\in{\mathbb{R}}^{(n+2)\times(n+1)}$. When deriving spectral methods the operators $J_{x}$ and $J_{y}$ play an important role as they can be interpreted as operators for “multiplication-by-$x$” and “multiplication-by-$y$,” respectively. That is,

[TABLE]

In other words, if f(x,y)={\mathbf{P}}(x,y)^{\top}\mbox{\boldmathf\unboldmath}, then J_{x}^{\top}\mbox{\boldmathf\unboldmath} and J_{y}^{\top}\mbox{\boldmathf\unboldmath} give the coefficients of $xf(x,y)$ and $yf(x,y)$, respectively.

2.3 Recurrences and the Clenshaw algorithm

For univariate orthogonal polynomials, the three-term recurrence encoded by a Jacobi operator can be used to construct the polynomials themselves at a specified point via forward substitution. Clenshaw’s algorithm is a closely related concept that allows the evaluation of a finite series expansion of univariate orthogonal polynomials at a point [4]. While it is common to interpret the three-term recurrence/Clenshaw’s algorithm as recursions, we prefer to interpret them as forward/backward substitution on a lower/upper triangular system associated to the Jacobi operator as this point-of-view facilitates generalization to the bivariate setting.

Let $p_{0}(x),p_{1}(x),\ldots,$ be a sequence of univariate orthogonal polynomials such that $p_{0}(x)=1$, and suppose that we wish to evaluate $f(x)=\sum_{k=0}^{N}a_{k}p_{k}(x)$ at $x_{*}\in\mathbb{R}$. Since $p_{0}(x),p_{1}(x),\ldots,$ satisfy a three-term recurrence of the form $b_{k}p_{k+1}(x)=(x-a_{k})p_{k}(x)-c_{k-1}p_{k-1}(x)$ for $k\geq 1$ [23, Thm. 3.2.1], we find that

[TABLE]

where $b_{0}p_{1}(x)=(a_{0}-x)p_{0}(x)$ and $e_{0}={{\left({1,0,\ldots,0}\right)}^{\top}}$.

Forward substitution on the lower triangular linear system in Eq. 9 allows one to evaluate $p_{k}(x_{*})$ for $k\geq 0$ from which one could evaluate $f(x_{*})=\sum_{k=0}^{N}a_{k}p_{k}(x_{*})$. For stability purposes, the Clenshaw algorithm evaluates expansions more directly and can be written as

[TABLE]

Therefore, the Clenshaw algorithm is equivalent to solving the upper triangular linear system $(L_{N}(x_{*}))^{\top}\boldsymbol{v}=\boldsymbol{a}$, followed by returning the first entry of $\boldsymbol{v}$. Since $L_{N}(x_{*})$ only has three nonzero subdiagonals, the algorithm requires $\mathcal{O}(N)$ operations to evaluate $f(x_{*})=\sum_{k=0}^{N}a_{k}p_{k}(x_{*})$.

The bivariate case is more involved. Given $(x_{*},y_{*})\in\mathbb{R}^{2}$, we would like to evaluate $f(x,y)=\sum_{n=0}^{N}\sum_{k=0}^{n}a_{n,k}p_{n,k}(x,y)$ at $(x_{*},y_{*})$, where (without loss of generality) we assume that $p_{0,0}(x,y)=1$. Since there are three-term recurrence relations in both $x$ and $y$ (see Eq. 7) we find that

[TABLE]

where $I_{m}$ is the $m\times m$ identity matrix and ${\bf 0}_{m\times 1}$ is the zero vector of length $m$. Unlike the univariate case, the system is not lower triangular and so we cannot immediately invert this system via forward recurrence to find $\mathbf{P}(x_{*},y_{*})$.

A reformulation that allows for inversion is to multiply the system to reduce the blocks above the diagonal in Eq. 11 to the identity. First, note that the blocks

[TABLE]

have full column rank for $n\geq 0$ [7, Theorem 3.3.4]. Therefore, $B_{n}$ has a left-inverse $B_{n}^{+}$ for $n\geq 0$ such that $B_{n}^{+}B_{n}=I_{n+2}$. It follows that an equivalent evaluation scheme can be designed from

[TABLE]

Since $\tilde{L}_{N}(x_{*},y_{*})$ is lower triangular we can construct $\mathbf{P}(x_{*},y_{*})$ via forward substitution.

Furthermore, a natural bivariate analogue of Clenshaw’s algorithm follows from writing

[TABLE]

Thus $f(x_{*},y_{*})$ can be evaluated by solving an upper triangular linear system using back substitution.

If $B_{n}^{+}$ are dense matrices for $n\geq 0$, then forward recurrence and Clenshaw’s algorithm require $\mathcal{O}(N^{3})$ operations. However, in the special case of Jacobi polynomials on the triangle, the matrices involved are sparse (see Section 3) and the complexity can be reduced to $\mathcal{O}(N^{2})$ operations, which is optimal.

2.4 Multiplication operators

The relations in Eq. 8 show that $J_{x}^{\top}$ and $J_{y}^{\top}$ are operators that represent “multiplication-by-$x$” and “multiplication-by-$y$”, respectively, in the bivariate orthogonal polynomial basis. Here, we combine these operators together to construct multiplication matrices that represent multiplication by a degree $d$ polynomial expanded as $q(x,y)=\sum_{n=0}^{d}\sum_{k=0}^{n}q_{n,k}p_{n,k}(x,y)$.

Suppose we are given a function $f(x,y)=\sum_{n=0}^{N}\sum_{k=0}^{n}a_{n,k}p_{n,k}(x,y)$, and wish to find the expansion coefficients of $g(x,y)=q(x,y)f(x,y)$, where the degree of $f$ and $q$ can differ. Using $J_{x}^{\top}$ and $J_{y}^{\top}$, we find that

[TABLE]

where the definition of $q(J_{x}^{\top},J_{y}^{\top})$ is444The matrix $M_{q}$ is the same as that studied in the literature on bivariate functions of matrices. More precisely, $M_{q}$ is denoted by $M_{q}=q\!\left\{J_{x}^{\top},J_{y}\right\}\!(I)$ in [12], where $I$ is the identity matrix and the missing transpose on the second argument is a matter of convention (see [12, Def. 2.1]).

[TABLE]

Since $J_{x}$ and $J_{y}$ are block-tridiagonal and each matrix-matrix product increases the block-bandwidth by one, we see that $q(J_{x}^{\top},J_{y}^{\top})$ is also a block-banded with upper and lower block-bandwidth $d$.

The expression in Eq. 13 is not ideal for computations when $d$ is moderately large because of the inherent ill-conditioning in the monomial basis. It is often computationally beneficial to expand $q(x,y)$ in a bivariate orthogonal polynomial expansion and evaluate $q(J_{x}^{\top},J_{y}^{\top})$ using an operator-valued analogue of Clenshaw’s algorithm [22, 28].

The operator-valued analogue of Clenshaw’s algorithm for evaluating $q(J_{x}^{\top},J_{y}^{\top})$ is equivalent to the expression:

[TABLE]

where $\mathcal{I}$ is an infinite identity matrix, $e_{0}$ is the first canonical unit vector, and $\boldsymbol{q}$ is the vector of coefficients for $q(x,y)$ in the bivariate orthogonal polynomial expansion. Here we use the Kronecker product denoted $\otimes$.

In general, $J_{x}$ and $J_{y}$ have dense blocks so that the total number of nonzero entries in the principal $N\times N$ block matrix of $q(J_{x}^{\top},J_{y}^{\top})$ is $\mathcal{O}(N^{3})$, and the complexity of constructing $q(J_{x}^{\top},J_{y}^{\top})$ using the operator-valued Clenshaw’s algorithm is $\mathcal{O}(N^{4})$ (where the total number of unknowns is $\mathcal{O}(N^{2})$). In the case of Jacobi polynomials on the triangle, the blocks of $J_{x}$ and $J_{y}$ are tridiagonal and there are only $\mathcal{O}(N^{2})$ nonzero entries, which can be calculated in optimal complexity.

3 Computing with Jacobi polynomials on the triangle

We now specialize the algorithmic ideas in Section 2 to Jacobi polynomials on the triangle (see Eq. 2). Since these polynomials have additional structure, more efficient algorithms can be designed.

We denote the Jacobi polynomials on the triangle that are orthogonal with respect to $x^{a}y^{b}(1-x-y)^{c}$ with $a,b,c>-1$ by

[TABLE]

and note that series expansions in the $P_{n,k}^{(a,b,c)}$ basis can be expressed as

[TABLE]

where $\boldsymbol{f}$ is the vector of $P_{n,k}^{(a,b,c)}$ coefficients for $f$. These expansion coefficients can be efficiently computed from samples of $f$ by a fast, backward stable algorithm [19, 20, 21]. We further denote the shifted Jacobi polynomials on the unit interval as

[TABLE]

where the alternative ordering of $a$ and $b$ helps to build analogies with the triangle case.

3.1 Conversion operators

An important property of Jacobi polynomials on the interval is that they have banded conversion operators, which translate between coefficients from expansion in $\smash{P_{n}^{(a,b)}}$ to $\smash{P_{n}^{(a+1,b)}}$ or $\smash{P_{n}^{(a,b+1)}}$. In terms of converting expansions between bases, we can express such conversions as

[TABLE]

where $S_{(a,b)}^{(a+1,b)}$ and $S_{(a,b)}^{(a,b+1)}$ are upper bidiagonal operators, with rational entries as given in [13, (18.9.5)].

Jacobi polynomials on the triangle have a similar property: we can increment either $a$, $b$, or $c$ in the expansion by one:

[TABLE]

Each of these operators are sparse: they have block-bandwidths $(0,1)$ with diagonal blocks for $S_{(a,b,c)}^{(a+1,b,c)}$ and upper bidiagonal blocks for $S_{(a,b,c)}^{(a,b+1,c)}$ and $S_{(a,b,c)}^{(a,b,c+1)}$. The entries are rational, and can be determined in closed form by the recurrence relationships in Corollary A.3.

3.2 Constructing Jacobi operators

For Jacobi polynomials, the recurrence relationships that give rise to tridiagonal Jacobi operators, representing multiplication by $x$, are well-known. However, the Jacobi operators can alternatively be derived via lower bidiagonal lowering operators $L_{(a,b)}^{(a-1,b)}$ and $L_{(a,b)}^{(a,b-1)}$ [13, 18.9.6] that represent multiplication by $x$ and $1-x$:

[TABLE]

Similarly, $1-x$ is equivalent to $S_{(a,b-1)}^{(a,b)}L_{(a,b)}^{(a,b-1)}$. In other words, the Jacobi operator corresponding to multiplication by $x$ can be constructed via

[TABLE]

Note that the product of a lower bidiagonal operator $L_{(a,b)}^{(a-1,b)}$ and an upper bidiagonal operator $S_{(a-1,b)}^{(a,b)}$ is a tridiagonal operator, as expected.

To construct the Jacobi operators $J_{x}$ and $J_{y}$ for Jacobi polynomials on the triangle, we first note that there exists three lowering operators that satisfy:

[TABLE]

where $z:=1-x-y$. We will use other indices to indicate multiple lowering in a row, e.g.,

[TABLE]

corresponds to multiplication by $xyz$, where the choice for navigating the parameter tree is arbitrary.

We can construct the Jacobi operators from the lowering operators via

[TABLE]

Note that the entries of the lowering operators can be determined by the recurrences in Corollary A.4, and they are sparse. In particular, they have block-bandwidths $(1,0)$ and diagonal blocks for $L_{(a,b,c)}^{(a-1,b,c)}$ and lower-bidiagonal blocks for $L_{(a,b,c)}^{(a,b-1,c)}$ and $L_{(a,b,c)}^{(a,b,c-1)}$. This block structure ensures that $J_{x}$ is block-tridiagonal with diagonal block and that $J_{y}$ is block-tridiagonal with tridiagonal blocks. Finally, $J_{x}$ and $J_{y}$ commute because the $L$ and $S$ operators commute:

[TABLE]

3.3 Implementation of Clenshaw’s algorithm and multiplication operators

We now exploit the sparsity structure of $J_{x}$ and $J_{y}$ to get an $\mathcal{O}(N)$ complexity Clenshaw algorithm. In particular, using the notation of Section 2.3, since $B_{n}^{x}$ is diagonal and $B_{n}^{y}$ is tridiagonal we can construct a simple left-inverse $B_{n}^{+}$. That is, we have the following structure:

[TABLE]

Let $B_{1}=B_{n}^{x}[0\!:\!n,0\!:\!n]$ denote the first $(n+1)\times(n+1)$ sub-block of $B_{n}^{x}$, let $b_{2}=B_{n}^{y}[n,n+1]$, and let

[TABLE]

Then, the following matrix is a pseudo-inverse of $B_{n}$:

[TABLE]

Note that $B_{n}^{+}$ can be applied to a vector in $\mathcal{O}(n)$ operations. When incorporation into Clenshaw’s algorithm described in Section 2.3, this gives an optimal $\mathcal{O}(N^{2})$ algorithm for evaluating functions. Furthermore, when incorporated into the construction of the multiplication operators (see Section 2.4), we find that one can construct multiplication operators in $\mathcal{O}(N^{2})$ operations.

3.4 Differentiation

Jacobi polynomials on the interval have banded recurrence relationships for their derivatives by incrementing both of the parameters, that is, we can represent

[TABLE]

where $D_{(a,b)}^{(a+1,b+1)}$ is zero except for the first super-diagonal [13, 18.9.15]. They also have banded recurrence relationship for their weighted derivatives that decrement the parameters:

[TABLE]

where $W_{(a,b)}^{(a-1,b-1)}$ is zero except for the first sub-diagonal [13, 18.9.16].

These properties translate to partial derivatives of Jacobi polynomials on the triangle. That is, we have555We have similar relationships for ${\partial\over\partial z}:={\partial\over\partial x}-{\partial\over\partial y}$, but we omit these for brevity as they are not needed.

[TABLE]

where the entries are derived in Corollary A.1. Both $D_{x,(a,b,c)}^{(a+1,b,c+1)}$ and $D_{y,(a,b,c)}^{(a,b+1,c+1)}$ are sparse: they are block super-diagonal, and their blocks are upper bi-diagonal and super-diagonal, respectively. Similarly, for weighted differentiation we have

[TABLE]

where the entries are derived in Corollary A.2. Both $W_{x,(a,b,c)}^{(a-1,b,c-1)}$ and $W_{y,(a,b,c)}^{(a,b-1,c-1)}$ are also sparse matrices as they are block sub-diagonal, and their blocks are lower bi-diagonal and sub-diagonal, respectively.

Combining differentiation and conversion appropriately allows us to represent more complicated differential operators. For example, the Laplacian can be expressed as an operator that takes coefficients in an ${\bf P}^{(0,0,0)}$ expansion to coefficients in an ${\bf P}^{(2,2,2)}$ expansion as follows:

[TABLE]

A simple calculation determines that this is also a sparse operator with block-bandwidths $(2,4)$ and blocks with bandwidths $(0,4)$, see the left figure in Fig. 2.

Similarly, we can express the Laplacian as an operator from coefficients in an $xy(1-x-y){\bf P}^{(1,1,1)}(x,y)$ expansion to coefficients in an ${\bf P}^{(1,1,1)}(x,y)$ expansion by using weighted derivatives and lowering operators:

[TABLE]

This is a sparse operator with block-bandwidths $(1,2)$ and blocks with bandwidths $(2,2)$, see the middle figure in Fig. 2.

Finally, variable coefficients can be constructed by combining lowering, conversion, and Jacobi operators. For example, the variable Helmholtz operator $\Delta+v(x,y)$ can be represented as

[TABLE]

This still leads to a sparse discretisation, where the block bandwidths depend on the degree of $v$, see the right figure in Fig. 2 for an example with $v(x,y)+xy^{2}$.

4 Solving linear PDEs with zero Dirichlet conditions

We now use the systematic approach to constructing sparse operators to solve PDEs. We construct the operators using BlockBandedMatrices.jl [14], which enables fast multiplication of block-banded matrices with banded blocks by building on BLAS. We then convert the representation to a SuiteSparse compatible sparse matrix format, for matrix factorization and solves.

4.1 Zero Dirichlet conditions

To solve PDEs with vanishing Dirichlet conditions, we use the weighted basis

[TABLE]

For higher order equations like the Biharmonic equation we consider vanishing Dirichlet and Neumann conditions using the weighted basis

[TABLE]

4.1.1 Example 1: Poisson equation

Consider Poisson’s equation on a triangle with zero Dirichlet conditions, i.e.,

[TABLE]

We reduce this equation to a truncation of

[TABLE]

where the coefficients of f(x,y)={\bf P}^{(1,1,1)}(x,y)^{\top}\mbox{\boldmathf\unboldmath} are determined using [19] as implemented in [21]. In Fig. 3 (left), we depict the solution for a specific choice of $f(x,y)$. In Figure 3 (right), we show the construction time666Timings are performed on an iMac 2017 with 3.8 GHz Intel Core i5, using Julia v1.0 compiled with MKL BLAS. Note that the default OpenBLAS is slower for banded matrix operations. of the matrix (using BlockBandedMatrices.jl), execution time for an LU factorization, and thesolve time (using SuiteSparse via Julia’s SparseArrays.jl). We observe that the construction requires an optimal $\mathcal{O}(N^{2})$ operations while the factorization and solution time are observed to cost an almost-optimal $\mathcal{O}(N^{3})$ operations. The same complexities are observed for PDEs on rectangles using a Chebyshev-based spectral method [9].

In Fig. 4, we show the norms of each block of calculated coefficients of the approximation for four right-hand sides with $N=999$. Note that the rate of decay in the coefficients is a proxy for the rate of convergence of the computed solution. The behavior of the right-hand side at the corners has an impact on the convergence rate; in particular, if $f$ and its derivatives vanish at the corners then we observe faster convergence of the solution. The behavior at the origin is particularly important as the Laplacian of the basis $xy(1-x-y)P_{n,k}^{(1,1,1)}(x,y)$ always vanishes at the origin. While we only observe algebraic convergence for the first three examples (that is, we do not achieve spectral convergence as $N\rightarrow\infty$), the rate of convergence is fairly fast, achieving machine precision accuracy when $f(x,y)$ vanishes at the origin with around 10,000 unknowns. Furthermore, the last example shows spectral convergence for a Gaussian bump function, which up to machine precision vanishes to all orders at the corners. Finally, over-resolving the solution does not result in the error plateauxing at machine precision, which means our discretization slightly improves the regularity of the data, similar to the ODE case in [15].

4.1.2 Example 2: Variable coefficient Helmholtz equation with forcing terms

Now, consider a variable coefficient Helmholtz equation with zero Dirichlet conditions, i.e.,

[TABLE]

We first approximate $v(x,y)$ by a polynomial [21] and then use the operator-valued Clenshaw’s algorithm to construct $v(J_{x}^{\top},J_{y}^{\top})$. We obtain the following discretization:

[TABLE]

where $J_{x}^{\top}$ and $J_{y}^{\top}$ are the Jacobi operators for ${\bf P}^{(1,1,1)}$ and

[TABLE]

In Fig. 5 we depict the solution for $k=100$ and plot the timings for construction, factorization, and solution for $k$ between 100 and 300, using polynomials of degree $2k$. The build time depends only on the discretization size, so we observe an $\mathcal{O}(k^{2})$ cost.

4.1.3 Example 3: The biharmonic equation

The same technique for constructing a sparse representation of the Laplacian $\Delta$ translates to the Biharmonic operator $\Delta^{2}$, though now we must use a basis that satisfies both zero Dirichlet and Neumann conditions. We can represent the Laplacian as a map from coefficients in an $x^{2}y^{2}(1-x-y)^{2}{\bf P}^{(2,2,2)}(x,y)$ expansion to coefficients in an ${\bf P}^{(0,0,0)}$ expansion by using weighted differentiation and lowering operators:

[TABLE]

Hence, the biharmonic operator can be sparsely represented as a map from coefficients in an $x^{2}y^{2}(1-x-y)^{2}{\bf P}^{(2,2,2)}(x,y)$ expansion to coefficients in an ${\bf P}^{(2,2,2)}$ expansion. This is simply given by $\Delta_{(0,0,0)}^{(2,2,2)}\Delta_{W^{2}}$.

In Fig. 6 we depict a solution to the biharmonic equation and show that the build time grows linearly with respect to the number of degrees of freedom employed to discretize the solution.

5 Nonzero Dirichlet conditions

To handle general nonzero Dirichlet boundary conditions, we wish to construct restriction operators that are sparse operators. To facilitate this, we use a basis where most elements of the basis vanish at the boundary. We take the weighted basis $x^{a}y^{b}(1-x-y)^{c}{\bf P}^{(a,b,c)}$, where $a,b,c$ are integers, and augment it with additional polynomials so that the basis can represent all bivariate polynomials. This is essentially the same procedure as in [10], but we do it in a way that preserves the sparsity of the restriction operators. Appendix B gives the definition of $Q_{n,k}^{(a,b,c)}(x,y)$, where $a,b,c\in\left\{{0,1}\right\}$, which is the basis we use to construct sparse discretizations. Here, most of $Q_{n,k}^{(1,b,c)}(0,y)$, most of $Q_{n,k}^{(a,1,c)}(x,0)$, and most of $Q_{n,k}^{(a,b,1)}(x,1-x)$ vanish.

Remark 1**.**

*Formally, $Q_{n,k}^{(a,b,c)}(x,y)$ can be thought of as $P_{n,k}^{(-a,-b,-c)}(x,y)$, which is made precise in [29] during the construction of the polynomials $J_{n,k}^{(a,b,c)}(x,y)$. However, the construction in [29] is normalized in a way that leads to underflow in double precision computing and we find it simpler to define our own basis $Q_{n,k}^{(a,b,c)}(x,y)$ in an ad hoc way. *

5.1 Derivative and conversion operators

Partial derivatives and conversion operators \mbox{\boldmathQ\unboldmath}^{(a,b,c)}(x,y) are similar to those derived for \mbox{\boldmathP\unboldmath}^{(a,b,c)}(x,y). Using the formulas in Corollary C.1, we can construct conversion operators that convert from one-edge bases to \mbox{\boldmathP\unboldmath}^{(0,0,0)}:

[TABLE]

Note that each operator is block upper bi-diagonal, with diagonal or upper bi-diagonal blocks. Similarly, Corollary C.3 derives sparse conversion operators from two-edge bases to one edge bases, which we denote by $\tilde{S}_{(1,1,0)}^{(1,0,0)}$, $\tilde{S}_{(1,1,0)}^{(0,1,0)}$, $\tilde{S}_{(1,0,1)}^{(1,0,0)}$, etc. Finally, Corollary C.5 derives sparse conversion operators from the three-edge basis to any of the two-edge bases, which we denote by $\tilde{S}_{(1,1,1)}^{(1,1,0)}$, $\tilde{S}_{(1,1,1)}^{(1,0,1)}$, and $\tilde{S}_{(1,1,1)}^{(0,1,1)}$. We can clearly compose these operators together to convert from, say, \mbox{\boldmathQ\unboldmath}^{(1,1,1)} to \mbox{\boldmathP\unboldmath}^{(0,0,0)}. For this purpose, we can define

[TABLE]

There are always several paths through the parameter space to convert one basis into another; however, any path that is chosen results in the same final conversion operator.

The same principle is true for derivatives, though for our purposes it suffices to restrict our attention to derivatives of the two-edge bases. Corollary C.7 gives us the entries for sparse (block-superdiagonal with at most bidiagonal blocks) operators that satisfy:

[TABLE]

General partial derivative operators can be constructed by combining conversion and derivative operators. For example, we can successfully construct the Laplacian from \mbox{\boldmathQ\unboldmath}^{(1,1,1)} to \mbox{\boldmathP\unboldmath}^{(1,1,1)} as

[TABLE]

This is a sparse operator with block-bandwidths and sub-blockbandwidths equal to $(1,4)$.

5.2 Restriction operators

The definitions of $Q_{n,k}^{(1,0,0)}(x,y)$, $Q_{n,k}^{(0,1,0)}(x,y)$, and $Q_{n,k}^{(0,0,1)}(x,y)$ each has a simple restriction formula to one of the three edges of the triangle:

[TABLE]

In other words, the restriction operator from expansion in $Q_{n,k}^{(1,0,0)}(x,y)$ to Legendre expansion on the edge from $(0,0)$ to $(0,1)$ is a block-banded operator, where the blocks themselves are very sparse: each block has precisely one nonzero entry. Similarly, the other two bases give restriction operators to the other edges. We denote these restriction operators as $R_{x}$, $R_{y}$, and $R_{z}$, respectively. They are given by

[TABLE]

where {\mbox{\boldmathP\unboldmath}}(x):={\bf P}^{(0,0)}(x)^{\top} are the shifted Legendre polynomials.

For the full Dirichlet operator, we need to restrict to all three edges. Thus we can construct restriction operators from \mbox{\boldmathQ\unboldmath}^{(1,1,1)} to the boundary, where the boundary basis are piecewise mapped Legendre polynomials. This restriction operator can be calculated by combining conversion and the one-edge restrictions as follows:

[TABLE]

This operator is also sparse as each component is a product of sparse operators.

5.3 The $\tau$-method

An issue we must deal with is boundary data with discontinuities at the corners. Consider, for example, the Laplace equation with Dirchlet conditions:

[TABLE]

If the boundary data has discontinuities, that is, $f(0,0)\neq g(0,0)$, $f(0,1)\neq h(0,1)$, or $g(1,0)\neq h(1,0)$, then the solution itself will have an arg-like singularity: e.g., near the origin we have the local behaviour

[TABLE]

Representing $u(x,y)$ by polynomials is therefore limited as they impose continuity. Other PDEs like the Helmholtz equation have similar before when the boundary data has discontinuities.

To overcome this issue, we adapt the Lanczos $\tau$-method, see [17] for an overview. The Lanczos $\tau$-method is a device to produce invertible systems for polynomial spectral methods by augmenting the equations with polynomial correction terms, that also provide error control by measuring the magnitude of the correction term. In our context we use it to capture discontinuities by augmenting the boundary data with corrections of the form

[TABLE]

That is, we add constants $\tau_{1}$ and $\tau_{2}$ to our discretisation:

[TABLE]

Now in our examples below we actually have mathematically continuous boundary data, however, round-off errors mean our boundary data is slightly discontinuous. The $\tau$ correction terms give a way of capturing this discontinuity without destroying the regularity of $u$. When the solution is resolved the $\tau$ terms are therefore negligible, and we can use the approximation of $u$ on its own.

Note that it is possible to add additional $\tau$ correction terms to make the system invertible, but this is a more technical task and hence we prefer to use a QR decomposition to solve the resulting rectangular linear system in a least squares sense. This does incur a substantial penalty, as SuiteSparse’s QR decomposition is significantly slower than its LU decomposition.

5.3.1 Example 4: Laplace’s equation

Consider Laplace’s equation with prescribed Dirichlet data:

[TABLE]

Expanding $f(x)$, $g(y)$, and $h(x)$ in Legendre series leads to a system of equations satisfied by $u$. That is,

[TABLE]

where $u$ are the coefficients of $u(x,y)$ in the basis \mbox{\boldmathQ\unboldmath}^{(1,1,1)}, $f$ are the Legendre coefficients of $f(y)$, $g$ are the Legendre coefficients of $g(x)$, and $h$ are the Legendre coefficients of $h(x)$. We augment this system with $\tau$ corrections, which are ultimately ignored in the approximation of the solution.

In Fig. 7 we plot the calculated coefficients for $N=999$ for three choices of boundary data: ${\rm e}^{x}\cos y$, $x^{2}$, and $x^{3}(1-x)^{3}(1-y)^{3}$. The first two examples exhibit algebraic decay, with the rate of decay dictated by the number of derivatives matching at the corners. The last example has a smooth solution (${\rm e}^{x}\cos y$ is harmonic) and we see that the algorithm achieves super-algebraic convergence, and is stable for large $N$. We also note that evaluating the approximation is exact to within an accuracy of $3\times 10^{-16}$ compared to the exact solution at the arbitrary point $(x,y)=(0.1,0.2)$.

5.3.2 Example 5: Transport equation

Nothing in this framework depends on the PDE being elliptic. Here, we consider the transport equation given by

[TABLE]

Information travels at a rate and direction dictated by $c$, and depending on its value we need either one or two edges to uniquely determine the solution. If $0\leq c\leq 1$ the solution is uniquely determined from the boundary on the bottom, and hence we use the basis $Q^{(0,1,0)}$. If $c>1$ then information is coming in from the right, so we use the basis $Q^{(0,1,1)}$ on the bottom and hypotenuse edges. If $c<0$ then information comes in from the left and we use the basis $Q^{(1,1,0)}$ on the bottom and left edges. In Fig. 8 we depict the three solutions.

6 Systems of PDEs

Systems of PDEs can be handled in a straightforward way by concatenating their blocks. As an example, we can solve the Poisson equation with Neumann conditions by re-expressing the PDE as a first-order system: writing $v=u_{x}$ expressed in the basis \mbox{\boldmathQ\unboldmath}^{(1,0,1)}, and $w=u_{y}$ expressed in the basis \mbox{\boldmathQ\unboldmath}^{(0,1,1)}, the system becomes

[TABLE]

6.0.1 Example 6: Helmholtz equation in a polygon

Note that being able to handle systems of PDEs in this manner also allows us to solve on polygonal domains that are partitioned into triangular elements. For example, consider the Helmholtz equation

[TABLE]

on the polygonal domain with the vertices $(0,0)$, $(1,0)$, $(1,1)$, $(0,2)$, $(0,1)$, and $(-1,1.5)$. We can decompose this domain into 4 triangles and represent the solution as well as its first derivatives in orthogonal polynomial expansions. This leads to a system of $4\times 3=12$ PDEs. We then impose continuity of the value and the normal derivative across the interfaces of each element, exploiting the fact that the restriction operator maps to the same basis of Legendre polynomials. (The orientation may be different, but reversing orientation of Legendre expansions corresponds to multiplying by a diagonal matrix that swaps the signs of every other coefficients.) We show the success of this approach in Fig. 9 for $k=10$ (left) and $k=20$ (right).

The discretization of the PDE system is sparse, and the complexity of building the matrices is an optimal $\mathcal{O}(N^{2})$ using degree $N$ polynomials within each element.

7 Conclusions

We have shown that bivariate orthogonal polynomials can lead to sparse discretizations of general linear PDEs on triangles with Dirichlet and Neumann boundary conditions. Instead of quadrature, we use sparse recurrence relationships combined with specialized linear algebra routines, allowing optimal complexity for building the linear systems. Multiple triangles can be patched together to solve PDEs on polygonal domains.

Another extension is to tetrahedra in 3D and higher. We expect this to be straightforward because the definitions of orthogonal polynomials on higher dimensional simplices is very similar to the 2D case. In 3D, we can use the following polynomials:

[TABLE]

which are orthogonal with respect to $x^{a}y^{b}z^{c}(1-x-y-z)^{d}$ on the unit 3D simplex. The most time-consuming part of such an extension is deriving the recurrences relationships. Note that in 3D and higher the sparsity of our construction is useful even for small discretisation sizes, as a degree $N$ dense discretisation (e.g. arising from collocation) would require calculating $\mathcal{O}(N^{6})$ entries, where the proposed construction would require an optimal $\mathcal{O}(N^{3})$ operations.

We used direct solvers via SuiteSparse to solve the resulting discretizations, which is fairly efficient with even millions of unknowns. However, to push the methodology further we will need robust iterative methods and the development of preconditioners. It is not yet clear how to design preconditioners in this setting.

Acknowledgments

This work began when the second author visited the first author at The University of Sydney. We are grateful for the travel support provided by The University of Sydney. We also thank Andrew Horning and Nicolas Boulle for carefully reading the draft and improving the text.

Appendix A Recurrence relationships for Jacobi polynomials on the triangle

Here, we outline the recurrence relationships for $P_{n,k}^{(a,b,c)}(x,y)$ that we employ, which were previously derived in [16, 29]. We define $z:=1-x-y$ and ${\partial\over\partial z}:={\partial\over\partial y}-{\partial\over\partial x}$.

Corollary A.1**.**

[16, Corollary 1]** The following recurrence relations for the partial derivatives hold:

[TABLE]

Corollary A.2**.**

[16, Corollary 2]** The following recurrence relations for the weighted partial derivatives hold:

[TABLE]

Corollary A.3**.**

[16, Corollary 3]** The following recurrence relations for conversions hold:

[TABLE]

Corollary A.4**.**

[16, Corollary 4]** The following recurrence relations for lowering operators hold:

[TABLE]

[TABLE]

[TABLE]

Appendix B Dirichlet basis definitions

Here, we define a basis, denoted by $Q_{n,k}^{(a,b,c)}(x,y)$, that we employ to impose general Dirichlet and Neumann boundary conditions. We construct $Q_{n,k}^{(a,b,c)}(x,y)$ by augmenting the weighted basis

[TABLE]

so that it spans all the polynomials with $a,b,c\in\left\{{0,1}\right\}$. Depending on the choice of $a$, $b$, and $c$, we obtain sparse restriction operators to one, two, or three edges of the triangle. We refer to this basis as the Dirichlet basis for its usefulness in solving PDEs with Dirichlet and Neumann boundary conditions.

B.1 One-edge Dirichlet basis

Definition B.1**.**

The following polynomials vanish at $x=0$ apart from when $k=n$:

[TABLE]

The following polynomials vanish at $y=0$ apart from when $k=0$:

[TABLE]

The following polynomials vanish at $z=0$ (i.e., $y=1-x$) apart from when $k=0$:

[TABLE]

The ordering is chosen so that the conversion operators derived below are upper triangular. Each basis has a simple restriction formula to the corresponding edge.

Proposition B.2**.**

Restriction operator to $x=0$:

[TABLE]

Restriction operator to $y=0$:

[TABLE]

Restriction operator to $z=0$:

[TABLE]

B.2 Two-edge Dirichlet basis

To handle two edges, consider first $x=0$ and $y=0$. As before, we wish to construct a basis that adds in the missing polynomials to $xyP_{n,k}^{(1,1,0)}(x,y)$ in a way that the restriction operators have the necessary structure. To do this, we select polynomials so that we can construct the conversion operator to expansions in the basis \mbox{\boldmathQ\unboldmath}^{(1,0,0)} and use the restriction operators we already have (see Corollary C.1).

Definition B.3**.**

The following polynomials vanish at $x=0$ and $y=0$ apart from when $k=0,n$:

[TABLE]

The following polynomials vanish at $x=0$ and $z=0$ apart from when $k=0,n$:

[TABLE]

The following polynomials vanish at $y=0$ and $z=0$ apart from when $k=0,1$:

[TABLE]

B.3 Three-edge Dirichlet basis

We finally get to three edges. Again, we want to choose the extra polynomials so that we can easily convert to any two-edge cases. The following does the trick:

Definition B.4**.**

The following polynomials vanish at $x=0$, $y=0$, and $z=0$ apart from when $k=0,1$, and $n$:

[TABLE]

Appendix C Dirichlet basis recurrence relationships

The following allows us to construct sparse conversion operators from the one-edge Dirichlet basis to the standard Jacobi polynomials on the triangle:

Corollary C.1**.**

The following recurrence relationships hold:

[TABLE]

Proof C.2**.**

*These are either immediate from definitions or are obtained by rearranging recurrence relationships found in Corollary A.4. *

The two-edge Dirichlet basis satisfy several sparse recurrence relationships.

Corollary C.3**.**

The following recurrence relationships hold:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Proof C.4**.**

*These are either immediate from definitions or are obtained by rearranging recurrence relationships found in Corollary A.4. *

The three-edge Dirichlet basis also satisfy several sparse recurrence relationships.

Corollary C.5**.**

The following recurrence relationships hold:

[TABLE]

[TABLE]

Proof C.6**.**

*These are either immediate from definitions or are obtained by rearranging recurrence relationships found in Corollary A.4. *

C.1 Recurrence relationships for the partial derivatives of the Dirichlet basis

We now turn to recurrence relationships for the partial derivatives of the Dirichlet basis, which are needed when imposing Neumann boundary conditions.

Corollary C.7**.**

The following recurrence relationships hold:

[TABLE]

Proof C.8**.**

The first three relations follow from the weighted partial differentiation relationships (see Corollary A.2). The fourth relation requires the additional property that

[TABLE]

which follows from [13, 15.5.6]. The fifth relationship also follows from Corollary A.2. For the sixth equation, if we define $t=y/(1-x)$, then the relation reduces to

[TABLE]

*and this expression follows from ${\cal L}_{2}^{\dagger}$ in [16, Lem. 1]. The last three relations follow from the same manipulation. *