A New Collocation Scheme Using Non-polynomial Basis Functions

Abstract

In this paper, we construct a set of non-polynomial basis functions from a generalised Birkhoff interpolation problem involving the operator: ${\mathscr {L}}_\lambda ={d^2}/{dx^2}-\lambda ^2$ with constant $\lambda .$ With a direct inverting the operator, the basis can be pre-computed in a fast and stable manner. This leads to new collocation schemes for general second-order boundary value problems with (i) the matrix corresponding to the operator ${\mathscr {L}}_\lambda$ being identity; (ii) well-conditioned linear systems and (iii) exact imposition of various boundary conditions. This also provides efficient solvers for time-dependent nonlinear problems. Moreover, we can show that the new basis has the approximability to general functions in Sobolev spaces as good as orthogonal polynomials.

Appendices

Formulas for $C_1$ and $C_2$ in Proposition 2.1

Corresponding to the typical boundary conditions in (2.3), we have the following formulas for the constant $C_1$ and $C_2$ in Proposition 2.1:

• if ${\mathscr {B}}_\pm [u]=u(\pm 1)=g_\pm ,$ then

  $$\begin{aligned} \begin{aligned}&C_1=\frac{1}{2\sinh (2\lambda )}\big \{g_++ {\mathcal {I}}_\lambda ^- [f](1)\big \} -\frac{e^{-2\lambda }}{2\sinh (2\lambda )}\big \{g_-+ {\mathcal {I}}_\lambda ^+ [f](-1)\big \},\\&C_2=-\frac{1}{2\sinh (2\lambda )}\big \{g_++ {\mathcal {I}}_\lambda ^- [f](1)\big \}+\frac{e^{2\lambda }}{2\sinh (2\lambda )}\big \{g_-+ {\mathcal {I}}_\lambda ^+ [f](-1)\big \}\big \}; \end{aligned} \end{aligned}$$

  (6.1)

• if ${\mathscr {B}}_\pm [u]=u'(\pm 1)=g_\pm ,$ then

  $$\begin{aligned} \begin{aligned}&C_1=\frac{\lambda ^{-1}}{2\sinh (2\lambda )}\big \{g_+-\lambda {\mathcal {I}}_\lambda ^- [f](1)\big \} -\frac{\lambda ^{-1}e^{-2\lambda }}{2\sinh (2\lambda )}\big \{g_-+\lambda {\mathcal {I}}_\lambda ^+ [f](-1)\big \},\\&C_2=\frac{\lambda ^{-1}}{2\sinh (2\lambda )}\big \{g_+-\lambda {\mathcal {I}}_\lambda ^- [f](1)\big \} -\frac{\lambda ^{-1}e^{2\lambda }}{2\sinh (2\lambda )}\big \{g_-+\lambda {\mathcal {I}}_\lambda ^+ [f](-1)\big \};\\ \end{aligned} \end{aligned}$$

  (6.2)

• if ${\mathscr {B}}_- [u]=u(-1)=g_-$ and ${\mathscr {B}}_+ [u]=u'(1)+ \eta \, u(1)=g_+,$ then

  $$\begin{aligned} \begin{aligned}&C_1=\frac{ g_++(\eta -\lambda ) {\mathcal {I}}_\lambda ^- [f](1)}{2\big (\lambda \cosh (2\lambda ) +\eta \sinh (2\lambda )\big )}- \frac{(\eta -\lambda ) e^{-2\lambda }\big \{g_- +{\mathcal {I}}_\lambda ^+ [f](-1)\big \}}{2\big (\lambda \cosh (2\lambda ) +\eta \sinh (2\lambda )\big )},\\&C_2=-\frac{g_++(\eta -\lambda ) {\mathcal {I}}_\lambda ^- [f](1)}{2\big (\lambda \cosh (2\lambda ) +\eta \sinh (2\lambda )\big )}+ \frac{(\eta +\lambda ) e^{2\lambda }\big \{g_- +{\mathcal {I}}_\lambda ^+ [f](-1)\big \}}{2\big (\lambda \cosh (2\lambda ) +\eta \sinh (2\lambda )\big )}. \end{aligned} \end{aligned}$$

  (6.3)

Jacobi Polynomials and Jacobi–Gauss–Lobatto Quadrature

Let $P_n^{(\alpha ,\beta )}(x)$ ($x\in [-1,1]$ and $\alpha ,\beta >-1$) be the Jacobi polynomial of degree n,  as normalized in [30]. We also refer to [30] for the following basic properties.

The Jacobi polynomials are eigenfunctions of the Sturm-Liouville equation

$$\begin{aligned} \begin{aligned} (x^2-1)\partial _x^2 P_n^{(\alpha ,\beta )}(x)+\big \{\alpha -\beta +(\alpha +\beta +2) x\big \}\partial _x P_n^{(\alpha ,\beta )}(x)=\lambda _n^{(\alpha ,\beta )} P_n^{(\alpha ,\beta )}(x), \end{aligned} \end{aligned}$$

(7.1)

where the corresponding eigenvalues are

$$\begin{aligned} \lambda _n^{(\alpha ,\beta )}=n(n+\alpha +\beta +1). \end{aligned}$$

(7.2)

The Jacobi polynomials are orthogonal with respect to the Jacobi weight function: $\omega ^{(\alpha ,\beta )}(x) = (1-x)^{\alpha }(1+x)^{\beta },$ namely,

$$\begin{aligned} \int _{-1}^1 {P}_n^{(\alpha ,\beta )}(x) {P}_{n'}^{(\alpha ,\beta )}(x) \omega ^{(\alpha ,\beta )}(x) \, \mathrm{d}x= \gamma _n^{(\alpha ,\beta )} \delta _{nn'}, \end{aligned}$$

(7.3)

where $\delta _{nn'}$ is the Dirac Delta symbol, and the normalization constant is given by

$$\begin{aligned} \gamma _n^{(\alpha ,\beta )} =\frac{2^{\alpha +\beta +1}\Gamma (n+\alpha +1)\Gamma (n+\beta +1)}{(2n+\alpha +\beta +1) n!\,\Gamma (n+\alpha +\beta +1)}. \end{aligned}$$

(7.4)

We have

$$\begin{aligned} P_n^{(\alpha ,\beta )}(x)=(-1)^n P_n^{(\beta ,\alpha )}(-x); \quad P_n^{(\alpha ,\beta )}(1)=\frac{\Gamma (n+\alpha +1)}{n!\,\Gamma (\alpha +1)}. \end{aligned}$$

(7.5)

Moreover, there holds the important derivative formula:

$$\begin{aligned} \partial _xP_n^{(\alpha ,\beta )}(x)=\frac{1}{2}(n+\alpha +\beta +1)P_{n-1}^{(\alpha +1,\beta +1)}(x),\;\;\; n\ge 1. \end{aligned}$$

(7.6)

We also use the following recurrent relation:

$$\begin{aligned} P_n^{(\alpha ,\beta )}(x)=a_n \partial _x P_{n-1}^{(\alpha ,\beta )}(x) +b_n \partial _x P_{n}^{(\alpha ,\beta )}(x)+c_n \partial _x P_{n+1}^{(\alpha ,\beta )}(x), \end{aligned}$$

(7.7)

where $a_1:=a_1^{(\alpha ,\beta )}=0$, and

$$\begin{aligned} a_n:=a_n^{(\alpha ,\beta )}=-\frac{2(n+\alpha )(n+\beta )}{(n+\alpha +\beta )(2n+\alpha +\beta )(2n+\alpha +\beta +1)},\quad n>1, \end{aligned}$$

(7.8a)

$$\begin{aligned} b_n:= b_n^{(\alpha ,\beta )}=\frac{2(\alpha -\beta )}{(2n+\alpha +\beta )(2n+\alpha +\beta +2)},\quad n\ge 1, \end{aligned}$$

(7.8b)

$$\begin{aligned} c_n:=c_n^{(\alpha ,\beta )}=\frac{2(n+\alpha +\beta +1)}{(2n+\alpha +\beta +1)(2n+\alpha +\beta +2)},\quad n\ge 1. \end{aligned}$$

(7.8c)

The Jacobi–Gauss–Lobatto (JGL) points $\big \{x_j=\xi _{N,j}^{(\alpha ,\beta )}\big \}_{j=0}^N$ (with $x_0=-1, x_N=1)$ are zeros of $(1-x^2) \partial _x P_{N}^{(\alpha ,\beta )}(x).$ Let $\big \{\omega _j=\omega _{N,j}^{(\alpha ,\beta )}\big \}_{j=0}^N$ be the corresponding JGL quadrature weights (cf. [27, Theorem. 3.27]). Then we have

$$\begin{aligned} \int _{-1}^1 \phi (x)\psi (x)\omega ^{(\alpha ,\beta )}(x)\,\mathrm{d}x=\sum _{j=0}^N \phi (x_j)\psi (x_j)\omega _j,\quad \forall \phi \cdot \psi \in {\mathcal {P}}_{2N-1}. \end{aligned}$$

(7.9)

Let $\big \{h_j:=h_{N,j}^{(\alpha ,\beta )}\big \}$ be the Lagrange interpolating basis polynomials associated with $\{x_j\}_{j=0}^{N},$ such that $h_j\in {\mathcal {P}}_{N}$ and $h_j(x_i)=\delta _{ij}.$ We have the representation

$$\begin{aligned} h_j(x)=\sum _{n=0}^N \frac{\omega _j}{\tilde{\gamma }_n} P_n^{(\alpha ,\beta )}(x_j) P_n^{(\alpha ,\beta )}(x), \end{aligned}$$

(7.10)

where

$$\begin{aligned} \tilde{\gamma }_n =\gamma _n^{(\alpha ,\beta )}, \;\;\; 0\le n\le N-1; \;\;\;\; \tilde{\gamma }_N=\Big ( 2+\frac{\alpha +\beta +1}{N} \Big )\gamma _N^{(\alpha ,\beta )}. \end{aligned}$$

Let ${\mathbb {I}}_N$ be the corresponding Lagrange interpolation operator, namely, ${\mathbb {I}}_N: C(\bar{\Lambda })\rightarrow {\mathcal {P}}_N$ such that for any $u\in C(\bar{\Lambda }),$

$$\begin{aligned} {\mathbb {I}}_N u(x_j)=u(x_j),\quad 0\le j\le N. \end{aligned}$$

(7.11)

Proof of Proposition 3.3

By the orthogonality (7.3) and (7.9),

$$\begin{aligned} \begin{aligned} \mu _{j}^n&=\frac{1}{\gamma _n^{(\alpha ,\beta )}}\int _{-1}^1 l_j(x) P_n^{(\alpha ,\beta )}(x)\,\omega ^{(\alpha ,\beta )}(x)\,\mathrm{d}x \\&=\frac{1}{\gamma _n^{(\alpha ,\beta )}}\Big \{l_j(-1)P_n^{(\alpha ,\beta )}(-1)\omega _0+ P_{n}^{(\alpha ,\beta )}(x_j)\omega _j + l_j(1) P_n^{(\alpha ,\beta )}(1)\omega _N \Big \}, \end{aligned} \end{aligned}$$

(8.1)

where we used the property: $l_j(x_i)=\delta _{ij}$ for $1\le i,j\le N-1$ (cf. (2.7)). Thus, it remains to derive the explicit formulas for $l_j(\pm 1).$ As the interior JGL points $\{x_j\}_{j=1}^{N-1}$ are zeros of $\partial _xP_N^{(\alpha ,\beta )}(x),$ we have

$$\begin{aligned} l_j(\pm 1)=\frac{\partial _xP_N^{(\alpha ,\beta )}(x)}{(x-x_j)\, \partial _x^2P_N^{(\alpha ,\beta )}(x_j)}\Bigg |_{x=\pm 1}. \end{aligned}$$

(8.2)

By (7.1),

$$\begin{aligned} \begin{aligned}&2(\beta +1) \partial _x P_N^{(\alpha ,\beta )}(-1)=-\lambda ^{(\alpha ,\beta )}_N P_N^{(\alpha ,\beta )}(-1), \;\; 2(\alpha +1) \partial _x P_N^{(\alpha ,\beta )}(1)=\lambda ^{(\alpha ,\beta )}_N P_N^{(\alpha ,\beta )}(1), \\&-(1-x^2_j)\partial _x^2 P^{(\alpha ,\beta )}_{N}(x_j)=\lambda ^{(\alpha ,\beta )}_N P^{(\alpha ,\beta )}_{N}(x_j),\quad 1\le j\le N-1. \end{aligned} \end{aligned}$$

(8.3)

A direct calculation from (8.2) and (8.3) leads to

$$\begin{aligned} l_j(-1)=-\frac{1-x_j}{2(\beta +1)}\frac{P_N^{(\alpha ,\beta )}(-1)}{P_N^{(\alpha ,\beta )}(x_j)},\quad l_j(1)=-\frac{1+x_j}{2(\alpha +1)}\frac{P_N^{(\alpha ,\beta )}(1)}{P_N^{(\alpha ,\beta )}(x_j)}. \end{aligned}$$

(8.4)

Thus, we obtain the desired formulas of $l_j(\pm 1)$ in (3.6).