1 Introduction↩︎

Finite element, finite difference, and spectral methods have long been employed in solving partial differential equations in one, two, and three space dimensions. Less has been done as the order of the equations increases. One exception to this is the application of Hermite splines of any order to solving two-point boundary problems of arbitrary even order [1]. The authors provide a priori error bounds on the finite element error but no computational results.

A dearth of models governed by higher-order equations may also explain the lack of interest. There are several exceptions to this. One involves blow-up solutions in one dimension for semilinear parabolic equations of order 2m+2 [2]–[4]. A second example comes from extensions of the seminal paper by Mindlin of his theory of deformation of elastic solids [5]. Both sixth- [6], [7] and eighth- [8] order equations have been derived from it to model the deformation of solids. A third example, for which fourth- and sixth-order equations have been studied, is a thin film model; see [9], [10] and references therein. Some preliminary analysis for a tenth-order extension of the equation has been made [11], [12].

Given the complications of increasing order it is appropriate to begin by examining equations in one space dimension. To that end consider two-point boundary value problems of the form \[\begin{align} f(x)&=(p_{m+1}(x)u^{(m+1)}(x))^{(m+1)}+\sum_{k=1}^m(p_k(x,u(x), \ldots, u^{(k-1)}(x))u^{(k)}(x))^{(k)}+p_0(x,u(x)), \, m \ge 0, \tag{1} \\ p_{m+1}(x)& > 0, \, x \in \Omega \equiv [a,b], \tag{2} \\ u_{a,k}&=u^{(k)}(a), \, u_{b,k}=u^{(k)}(b), \, k=0, 1, \ldots, m, \tag{3} \end{align}\] where the functions \(f\), \(p_k, \, k=0, 1, \ldots m+1\), and constants \(u_{a,k}, \, u_{b,k}, \, k=0, 1, \ldots m\), are such that solutions \(u(x) \in C^{\infty}(\Omega)\) exist and are unique. The weak form of 1 , with \(u(x) \in H^{(m+1)}(\Omega)\), is given by \[\begin{align} \int_a^b&f(x)v(x)dx=\int_a^b(-1)^{m+1}u^{(m+1)}(x)v^{(m+1)}(x)dx \nonumber \\ &+\int_{a}^b\Bigg(\sum_{k=1}^m(-1)^kp_k(x,u(x), \ldots, u^{(k-1)}(x))u^{(k)}(x)v^{(k)}(x) +p_0(x,u(x))v(x)\Bigg)dx, \, \forall v(x) \in H_0^{m+1}(\Omega). \end{align}\] The finite element solution is obtained by approximating \(H_0^{m+1}(\Omega)\) and \(H^{m+1}(\Omega)\) by finite dimensional subspaces. The choice of finite dimensional subspace determines the accuracy of the solution. The results shown extend to time-dependent problems and to equations on tensor-product domains in two and three space dimensions.

This approach has some advantages over the alternative that replaces the high-order equation with a coupled system of lower-order equations. If the goal is to control the error in \(H^{m+1}(\Omega)\) it is more natural to leave the equation in its original form. The same is true regarding the distribution of the additional degrees of freedom if a basis of order higher than \(2m+1\) is employed. For example, if the equation is replaced by a coupled system of second order equations and the basis order remains \(p > 2m+1\) across all equations then the system contains an additional \((p-2m-1)(m+1)\) unknowns per element.

One of the difficulties in solving 1 directly using finite elements is the determination of nodal basis functions to preserve \(C^m\) continuity across element boundaries. Linear and cubic Hermite splines are an obvious and common choice for \(m=0\) and \(m=1\), respectively. These splines have been extended from \(m=2\) to \(m=11\) [13]–[17] by deriving the coefficients from the continuity conditions. A general formula for any order can be derived from [18]. An alternative approach that involves one polynomial per interval based at its midpoint which matches solutions and \(m+1\) derivatives at the element boundaries is employed in [19] in a staggered grid solver for hyperbolic equations.

The complementary interior basis functions are also presented along with their key properties. These interior basis functions are formed from the product of lower-order ultraspherical polynomials multiplied by binomial powers of degree \(m+1\) associated with the element boundaries. They turn out to be natural extensions of Jacobi polynomials. Finally, it is shown that the two families of polynomials are linked through the derivatives of the Hermite splines.

The building blocks of Section 2 are used to construct the Hermite-Jacobi interpolant in Section 3. This is first carried out on a single interval \(\Delta=[M-h/2,M+h/2]\). Interpolation is done by matching the underlying function and its first \(m\) derivatives at \(M-h/2\) and \(M+h/2\) and at the scaled roots of the paired ultraspherical polynomial. A formula for the error in the interpolant is derived and extended to errors in \(L^2(\Delta)\) and \(H^k(\Delta), \, k=1, 2, \ldots m+1\) norms. In proving the error formula, generalizations of the Lagrange interpolants are introduced. The Hermite-Jacobi interpolant is extended to a grid that discretizes the interval \(\Omega=[a,b]\) along with the \(L^2(\Omega)\) and \(H^k(\Omega)\) errors. Superconvergence and the asymptotic equivalence of the Hermite-Jacobi interpolant are proved for a linear equation in Section 4. Computations for a model problem demonstrate the theory in Section 5 with conclusions and next steps presented in Section 6.

2 Arbitrary-Order Hermite and Jacobi Polynomials: Their Basis Properties and Interrelationships↩︎

Consider the two families of polynomials on the two canonical intervals \(\Omega_0=[0,1]\) and \(\Omega_1=[-1,1]\), respectively. On \(\Omega_0\) consider the pair of polynomials of order \(2m+1\), \[\begin{align} {\cal H}_l^m(\xi)&=\frac{(1-\xi)^{m+1}\xi^l}{l!} \sum_{n=0}^{m-l} S_n^m \xi^n, \, \tag{4} \\ {\cal I}_l^m(\xi)&=(-1)^l{\cal H}_l^m(1-\xi), \, l=0, 1, \ldots m, \, \xi \in \Omega_0, \tag{5} \end{align}\] where \[S_n^m=\frac{(m+n)!}{m!n!} \label{eq:defsmn}.\tag{6}\] These pairs of polynomials for each \(m\) serve as the left and right nodal basis functions, respectively. As will be shown below they ensure \(C^{m+1}\) continuity across element boundaries. On \(\Omega_1\) consider the polynomials of order \(p\) with zeros of order \(m+1\) at \(-1\) and \(1\), \[{\cal J}_p^m(\mu)=J_{m,p}(\mu^2-1)^{m+1}P_{p-m-1}^{(m+1)}(\mu)\equiv(\mu^2-1)^{m+1}{\cal J}_{p-2m-2}^{-m-2}(\mu), \, p \ge 2m+2, \, m \ge -1, \, \mu \in \Omega_1, \label{eq:ultrasphericaldef}\tag{7}\] where \(P_{p-m-1}(\mu)\) is the Legendre polynomial of degree \(p-m-1\) and \[J_{m,p}=\frac{2^{p-m-2}(p-m-2)!(p-2m-2)!}{(2p-2m-3)!}, \, p \ge 2m+2, \label{eq:h2}\tag{8}\] is chosen so that \({\cal J}_p^m(\mu)\) is monic. It also follows that \[{\cal J}_{p-2m-2}^{-m-2}(\mu) \equiv J_{m,p}P_{p-m-1}^{(m+1)}(\mu). \label{eq:JNormalize}\tag{9}\] The polynomials \({\cal J}_p^m(\mu)\), \(m \ge 0\) serve as the interior basis functions on an element. As will be shown, their counterparts, \({\cal J}_{p-2m-2}^{-m-2}(\mu)\) turn out to be the ultraspherical polynomials \({\cal C}_{p-2m-2}^{(\alpha)}(\mu)\) with \(2\alpha=2m+3\) in monic form. While both \({\cal J}_p^m(\mu)\) and \({\cal J}_{p-2m-2}^{-m-2}(\mu)\) form orthogonal sets of polynomials, the latter is a complete set while the former is not except when \(m=-1\) when they are equal and equal to Legendre polynomials in monic form.

The key properties of these two families are laid out in the next two theorems. Their relationship is presented in the third theorem.

The structure of these Hermite splines is best illustrated by considering two examples of the \({\cal H}_l^m(\xi)\). For the first example let \(m=1\). Then \[{\cal H}_0^1(\xi)=(1-\xi)^2(1+2\xi), \, {\cal H}_1^1(\xi)=\xi(1-\xi)^2, \nonumber\] are the traditional cubic Hermite basis pair on [0,1] at the left end. The second, with \(m=5\), produces \[\begin{align} {\cal H}_0^{5}(\xi)=(1-\xi)^6(252\xi^5+126\xi^4+56\xi^3+21\xi^2+6\xi+1), \nonumber \\ {\cal H}_1^{5}(\xi)=\xi(1-\xi)^6(126\xi^4+56\xi^3+21\xi^2+6\xi+1), \nonumber \\ {\cal H}_2^{5}(\xi)=\xi^2(1-\xi)^6(56\xi^3+21\xi^2+6\xi+1)/2!, \nonumber \\ {\cal H}_3^{5}(\xi)=\xi^3(1-\xi)^6(21\xi^2+6\xi+1)/3!, \nonumber \\ {\cal H}_4^{5}(\xi)=\xi^4(1-\xi)^6(6\xi+1)/4!, \nonumber \\ {\cal H}_5^{5}(\xi)=\xi^5(1-\xi)^6/5!. \nonumber \end{align}\] As \(l\) increases one additional term is removed from the sum and the result is divided by \(l\), that is, \[{\cal H}_l^m(\xi)=({\cal H}_{l-1}^m(\xi)-S_{m-l+1}^m\xi^{m-l+1})\xi/l, \, l > 0. \label{eq:hermiterecursion}\tag{10}\] In [20] (see section 2.1.5, p. 53) there is a similar recursion formula but in reverse, that is, starting with \({\cal H}_m^m(\xi)\) but unlike 10 it involves the computation of a derivative of an auxiliary polynomial and thus is not explicit.

Having demonstrated the basis properties of the Hermite splines consider now the second family defined by 7 8 .

Two lemmas are needed before proceeding to the proof of the theorem.

Proof (of Theorem 2). Equation ?? follows directly by differentiating 8 . From 7 8 it follows that \[J_{m-1,p-1}P_{p-m-1}^{(m)}(\mu)={\cal J}_{p-2m-1}^{-m-1}(\mu). \label{eq:h2a}\tag{12}\] Differentiating 12 and using 7 and ?? leads to ?? . To establish ?? , equation 7 implies that \[\begin{align} ({\cal J}_{p}^{m}(\mu))^{\prime}&=J_{m,p}(\mu^2-1)^{m+1}P_{p-m-1}^{(m+2)}(\mu)+J_{m,p}2\mu(\mu^2-1)^{m}(m+1)P_{p-m-1}^{(m+1)}(\mu) \nonumber \\ &=-J_{m,p}(\mu^2-1)^{m}((1-\mu^2)P_{p-m-1}^{(m+2)}(\mu)-2\mu(m+1)P_{p-m-1}^{(m+1)}(\mu)) \nonumber \\ &=J_{m,p}(\mu^2-1)(p)(p-2m-1)P_{p-m-1}^{(m)}(\mu) \nonumber \\ &=p{\cal J}_{p-1}^{m-1}(\mu), \label{eq:jintprime} \end{align}\tag{13}\] where ?? and ?? have been used. Equation ?? follows by integrating 13 together with ?? . From [21] (see equation 5, section 87) \[(2p+1)P_p(\mu)=P_{p+1}^{\prime}(\mu)-P_{p-1}^{\prime}(\mu). \label{eq:legderiva}\tag{14}\] Integrating 14 over \([-1,\mu]\) and using 7 and ?? results in \[\int_{-1}^{\mu}{\cal J}_p^{-1}(\rho)d\rho=\frac{1}{p+1}{\cal J}_{p+1}^{-1}(\mu)-\frac{p}{(2p+1)(2p-1)}{\cal J}_{p-1}^{-1}(\mu),\] or \[{\cal J}_{p+1}^0(\mu)={\cal J}_{p+1}^{-1}(\mu)-\frac{p(p+1)}{(2p+1)(2p-1)}{\cal J}_{p-1}^{-1}(\mu). \label{eq:legderivb}\tag{15}\] Replacing \(p+1\) with \(p\) in 15 implies that ?? holds when \(m=0\). By induction assume ?? is true for \(m=0, 1, \ldots, M\). Integrating ?? with \(m=M\) yields \[\int_{-1}^{\mu}{\cal J}_p^{\bar M}(\rho)d\rho=\int_{-1}^{\mu}{\cal J}_{p}^{{\bar M}-1}(\rho)d\rho-\frac{p(p-1)}{(2p-2{\bar M}-1)(2p-2{\bar M}-3)}\int_{-1}^{\mu}{\cal J}_{p-2}^{{\bar M}-1}(\rho)d\rho, \, p \ge 2{\bar M}+2,\] or using ?? \[\frac{1}{p+1}{\cal J}_{p+1}^{{\bar M}+1}(\mu)=\frac{1}{p+1}{\cal J}_{p+1}^{\bar M}(\mu)-\frac{p}{(2p-2{\bar M}-1)(2p-2{\bar M}-3)}{\cal J}_{p-1}^{\bar M}(\mu), \, p \ge 2{\bar M}+3. \label{eq:legderivc}\tag{16}\] Replacing \(p+1\) by \(p\) in 16 yields ?? with \(m={\bar M}+1\). If \({\bar C}_p^{(\alpha)}(\mu)\) is the monic form of the ultraspherical polynomial then \[{\bar {\cal C}}_p^{(\alpha)}(\mu)=\frac{(2\alpha-1)!p!2^{p-1}((2p+2\alpha-3)/2)!}{(2\alpha+2p-2)!((2\alpha-1)/2)!}{\cal C}_p^{(\alpha)}(\mu). \label{eq:monicultra}\tag{17}\] follows from [22] (see equation 22.3.4 p. 775) along with the three-term recurrence relation \[\begin{align} {\bar {\cal C}}_{p+1}^{(\alpha)}(\mu)&=\mu{\bar {\cal C}}_p^{(\alpha)}(\mu)-\frac{p(p+2m+2)}{(2p+2m+3)(2p+2m+1)}{\bar {\cal C}}_{p-1}^{(\alpha)}(\mu), \, p \ge 1, \tag{18} \\ {\bar {\cal C}}_0^{(\alpha)}&=1, \, {\bar {\cal C}}_1^{(\alpha)}=\mu, \tag{19} \end{align}\] [22] (see equation 22.7.3 p. 782). Differentiating the three-term recurrence formula for Legendre polynomials in [21] (see equation 2 section 89) yields \[pP_p^{\prime}(\mu)=(2p-1)P_{p-1}(\mu)+(2p-1)\mu P_{p-1}^{\prime}(\mu)-(p-1)P_{p-2}^{\prime}(\mu). \label{eq:legendrerecurrence}\tag{20}\] Applying 9 and ?? with \(m=0\) to 20 leads to \[\frac{p}{J_{0,p+1}}{\cal J}_{p-1}^{-2}(\mu)=(2p-1)P_{p-1}(\mu)+\mu \frac{2p-1}{J_{0,p}}{\cal J}_{p-2}^{-2}(\mu)-\frac{p-1}{J_{0,p-1}}{\cal J}_{p-3}^{-2}(\mu). \label{eq:jreducerecurrence}\tag{21}\] From [21] (equation 6 section 87), ?? , and 9 it is also true that \[P_{p-1}(\mu)=\frac{1}{p}\left( \frac{1}{J_{0,p+1}}{\cal J}_{p-1}^{-2}(\mu)-\mu\frac{1}{J_{0,p}}{\cal J}_{p-2}^{-2}(\mu)\right). \label{eq:diffrecurrence}\tag{22}\] Substituting 22 into 21 and using 8 gives \[{\cal J}_{p+1}^{-2}(\mu)=\mu{\cal J}_{p}^{-2}(\mu)-\frac{p(p+2)}{(2p+3)(2p+1)}{\cal J}_{p-1}^{-2}(\mu),\] which is the same as 18 with \(m=0\) or \(\alpha=3/2\). Since 8 and 9 imply that \({\cal J}_0^{-2}(\mu)=1\) and \({\cal J}_1^{-2}(\mu)=\mu\) it follows that ?? holds for \(m=0\). From [23] \[{\bar {\cal C}}_p^{(\alpha) \prime}(\mu)=p{\bar {\cal C}}_{p-1}^{(\alpha+1)}. \label{eq:ultramonicderiv}\tag{23}\] Then ?? , ?? with \(m=0\), and 23 imply that ?? holds for \(m > 0\) by successive differentiation. This implies that ?? holds for \(m \le -1\). To complete the proof of ?? when \(m > -1\) first note that from 7 that the first two interior basis functions are \[{\cal J}_{2m+2}^m(\mu)=(\mu^2-1)^{m+1}, \, {\cal J}_{2m+3}^m(\mu)=\mu(\mu^2-1)^{m+1}. \label{eq:firsttwointerior}\tag{24}\] Integrating both sides of the three-term recurrence formula for Legendre polynomials in [21] (see equation 2, section 89), applying integration parts to the first term on the right, and using ?? and ?? yields \[\begin{align} &p\int_{-1}^{\mu}P_p(\rho)d\rho=\frac{(2p-1)(2p-3)!}{2^{p-2}p!(p-2)!}\mu {\cal J}_p^0(\mu)-\frac{(2p-1)(2p-3)!}{2^{p-2}p!(p-2)!}\int_{-1}^{\mu}{\cal J}_p^0(\rho)d\rho \nonumber \\ &-(p-1)\int_{-1}^{\mu}P_{p-2}(\rho)d\rho, \, p \ge 2. \label{eq:phi1b} \end{align}\tag{25}\] Applying ?? with \(m=0\) to the second integral in 25 results in \[\begin{align} &p\int_{-1}^{\mu}P_p(\rho)d\rho=\frac{(2p-1)(2p-3)!}{2^{p-2}p!((p-2)!}\mu {\cal J}_p^0(\mu) \nonumber \\ &-\frac{(2p-1)(2p-3)!}{2^{p-2}p!((p-2)!}\int_{-1}^{\mu}\left( {\cal J}_p^{-1}(\rho)-\frac{p(p-1)}{(2p-1)(2p-3)}{\cal J}_{p-2}^{-1}(\rho)\right)d\rho-(p-1)\int_{-1}^{\mu}P_{p-2}(\rho)d\rho, \, p \ge 2. \label{eq:phi1f} \end{align}\tag{26}\] Combining like terms in 26 and applying ?? to the first and last terms together with 7 8 and \(m=-1\) yields \[\frac{(2p-3)!(2p-1)}{2^{p-2}p!(p-2)!}{\cal J}_{p+1}^0(\mu)=\frac{(2p-1)(2p-3)!}{2^{p-2}p!((p-2)!}\mu{\cal J}_p^0(\mu)-\frac{(p-2)(2p-5)!}{2^{p-3}(p-1)!(p-3)!}{\cal J}_{p-1}^0(\mu), \, p \ge 3, \label{eq:phi1d}\tag{27}\] where \(p\) is now greater than 3 since \({\cal J}_2^0\) is the lowest degree interior polynomial that satisfies ?? . Dividing both sides of 27 by the coefficient of the first term gives \[{\cal J}_{p+1}^0(\mu)=\mu{\cal J}_p^0(\mu)-\frac{2^{p-2}p!(p-2)!(p-2)(2p-5)!}{(2p-3)!(2p-1)2^{p-3}(p-1)!(p-3)!}{\cal J}_{p-1}^0(\mu), \, p \ge 3,\] which simplifies to \[{\cal J}_{p+1}^0(\mu)=\mu{\cal J}_p^0(\mu)-\frac{p(p-2)}{(2p-1)(2p-3)}{\cal J}_{p-1}^0(\mu), \, p \ge 3. \label{eq:phi1e}\tag{28}\] Equation 28 is true for \(p=2\) since in that case the second term on the right-hand side is zero and since 24 holds. When \(p=1\), 28 , using 24 can be written as \[\mu{\cal J}_1^0(\mu)-{\cal J}_0^0(\mu)={\cal J}_2^0(\mu)\equiv(\mu^2-1),\] which is satisfied if \({\cal J}_1^0(\mu)=\mu\) and \({\cal J}_0^0(\mu)=1\). Thus, ?? holds for \(m=0\). By induction assume true for \(m=0, 1, \ldots, {\hat{M}}\). Integrating ?? and using integration by parts on the second term together with ?? results in \[\begin{align} &\frac{1}{p+1}{\cal J}_{p+1}^{{\hat{M}}+1}(\mu)=\frac{1}{p}\mu{\cal J}_p^{{\hat{M}}+1}(\mu)-\frac{1}{p}\int_{-1}^{\mu}{\cal J}_p^{{\hat{M}}+1}(\rho)d\rho -\frac{(p-2{\bar M}-3)}{(2p-2{\bar M}-3)(2p-2{\bar M}-5)}{\cal J}_{p-1}^{{\hat{M}}+1}(\mu), \nonumber \\ &p \ge 2{\bar M}+5. \label{eq:phi1g} \end{align}\tag{29}\] Substituting from ?? for the integral in 29 , combining like terms, and finally applying ?? to the first and last terms yields \[\left(1+\frac{1}{p} \right) \frac{1}{p+1}{\cal J}_{p+1}^{{\hat{M}}+1}(\mu)=\frac{1}{p}\mu{\cal J}_p^{{\hat{M}}+1}(\mu)-\frac{(p-2{\bar M}-4)}{(2p-2{\bar M}-3)(2p-2{\bar M}-5)}{\cal J}_{p-1}^{{\hat{M}}+1}(\mu),\] or \[{\cal J}_{p+1}^{{\hat{M}}+1}(\mu)=\mu{\cal J}_p^{{\hat{M}}+1}(\mu)-\frac{p(p-2{\bar M}-4)}{(2p-2{\bar M}-3)(2p-2{\bar M}-5)}{\cal J}_{p-1}^{{\hat{M}}+1}(\mu), \, p \ge 2{\bar M}+5, \label{eq:phi1h}\tag{30}\] which confirms the three-term recurrence relation for \({\hat{M}}+1\). That 30 holds for \(p=2{\bar M}+4\) follows from 24 since the second term on the right-hand side is zero. To show that \({\cal J}_1^{{\bar M}+1}(\mu)=\mu\) and \({\cal J}_0^{{\bar M}+1}(\mu)=1\) rewrite 30 with \(p=2{\bar M}+3\) as \[\mu{\cal J}_{2{\bar M}+3}^{{\bar M}+1}(\mu)+\frac{1}{2{\bar M}+1}{\cal J}_{2{\bar M}+2}(\mu)={\cal J}_{2{\bar M}+4}^{{\bar M}+1}(\mu)\equiv(\mu^2-1)^{{\bar M}+2}, \label{eq:phi1i}\tag{31}\] where 24 has been used. Assume \[{\cal J}_p^{{\bar M}+1}(\mu)=\sum_{k=0}^p j_{p,k}^{{\bar M}+1}\mu^k, \, p \le 2{\bar M}+3. \label{eq:jacobiextn}\tag{32}\] Plugging 32 into 31 and equating like powers of \(\mu\) on both sides yields \[j_{p,p}^{{\bar M}+1}=1, \, j_{p,p-1}^{{\bar M}+1}=0. \label{eq:phi1j}\tag{33}\] Substituting \(p-1\) for \(p\) in 32 and using 30 along with 33 implies that \(j_{p-1,p-1}^{{\bar M}+1}=1\) and \(j_{p-1,p-2}^{{\bar M}+1}=0\). Induction on \(p\) implies that \({\cal J}_1^{{\bar M}+1}(\mu)=\mu\). Reversing directions then implies that \({\cal J}_p^{{\bar M}+1}(\mu)\) is odd or even as \(p\) is odd or even, respectively. From 31 it follows that \(j_{p-1,0}^{{\bar M}+1}=(2{\bar M}+1)(-1)^{{\bar M}+2}\). Continuing to move down in order and using the odd/even order of the \({\cal J}_p^{{\bar M}+1}(\mu)\) leads to \[j_{2,0}^{{\bar M}+1}=(-1)^{{\bar M}+2}\frac{1\cdot(-1)\cdot(-3) \cdots (-2{\bar M}+1)}{3\cdot5\cdots(2{\bar M}-1)(2{\bar M}+1)}=\frac{1}{2{\bar M}+1},\] which in turn implies that \({\cal J}_0^{{\bar M}+1}(\mu)=1\) using 30 with \(p=2\). Thus, ?? is established for all \(m\). To prove that the \({\cal J}_p^m(\mu)\) are orthogonal with respect to the weight function \((1-\mu^2)^{-(m+1)}\) from 7 it follows that \[\int_{-1}^1 \frac{{\cal J}_p^m(\rho){\cal J}_q^m(\rho)}{(1-\rho^2)^{m+1}}d\rho=\int_{-1}^1{\cal J}_{p-2m-2}^{-m-2}(\rho){\cal J}_{q-2m-2}^{-m-2}(\rho)(1-\rho^2)^{m+1}d\rho. \label{eq:jorthoa}\tag{34}\] Substituting into 34 from ?? yields \[\begin{align} \int_{-1}^1 \frac{{\cal J}_p^m(\rho){\cal J}_q^m(\rho)}{(1-\rho^2)^{m+1}}d\rho&=\frac{((2m+2)!)^2(p-2m-2)!(q-2m-2)!2^{p+q-4m-6}(p-m-2)!(q-m-2)!}{(2p-2m-3)!(2q-2m-3)!((m+1)!)^2} \nonumber \\ &\times \int_{-1}^1C_{p-2m-2}^{(\alpha)}C_{q-2m-2}^{(\alpha)}(1-\rho^2)^{\alpha-1/2}d\rho. \label{eq:orthob} \end{align}\tag{35}\] This establishes ?? in the case \(p \ne q\) from the orthogonality of the ultraspherical polynomials. Since from [22] (see equation 22.2.3 p. 774) \[\begin{align} \int_{-1}^1(C_{p-2m-2}^{(\alpha)}(\rho))^2(1-\rho^2)^{\alpha-1/2}d\rho&=\frac{\pi2^{1-2\alpha}\Gamma(p-2m-2+2\alpha)}{(p-2m-2)!(p-2m-2+\alpha)(\Gamma(\alpha))^2} \nonumber \\ &=\frac{2^{2m+1}p!(m!)^2}{(p-2m-2)!(2p-2m-1)((2m+1)!)^2}, \end{align}\] equation 35 simplifies to \[\begin{align} \int_{-1}^1 \frac{({\cal J}_p^m(\rho))^2}{(1-\rho^2)^{m+1}}d\rho&=\frac{2^{2m+1}p!(m!)^2}{(p-2m-2)!(2p-2m-1)((2m+1)!)^2} \nonumber \\ &\times \frac{((2m+2)!)^2((p-2m-2)!)^22^{2p-4m-6}((p-m-2)!)^2}{((2p-2m-3)!)^2(m+1)!}, \end{align}\] which further simplifies to ?? . 0◻

Remark: The polynomials \({\cal J}_p^m(\mu)\) are natural extensions of the Jacobi polynomials \(P_p^{(-m-1,-m-1)}(\mu)\), \(m \ge 0\), and are therefore also extensions of the ultraspherical polynomials \(C_p^{(\alpha)}, \, m \le -2,\) based on extending the three-term recurrence relation including the same two starting polynomials.

The Hermite splines and the polynomials \({\cal J}_p^m(\mu)\) defined by 7 8 are connected as demonstrated in the next theorem.

3 The Hermite-Jacobi Interpolant↩︎

The Hermite-Jacobi interpolant (the latter name is chosen since the \({\cal J}_p^m(\mu)\) are extensions of Jacobi polynomials) is constructed from the two families of polynomials introduced in the previous section. The Hermite splines, 4 6 , serve as the basis functions at the nodes, while the subfamily of the extended Jacobi polynomials, 7 8 , are the basis functions in the interior. Consider the interval \(\Delta\equiv[M-h/2,M+h/2]\) and let \(u(x) \in C^{\infty}({\bar \Delta})\). The approximation space \(S^{p,\Delta}\) consists of polynomials of degree \(p \ge 2m+1\). To that end define \[\begin{align} \Phi_{i,L,\Delta}^m(x)&={\cal H}_i^m((x-(M-h/2))/h), \, i=0, \cdots , m, \tag{36} \\ \Phi_{i,R,\Delta}^m(x)&={\cal I}_i^m((x-(M-h/2))/h), \, i=0, \cdots , m, \tag{37} \\ \Upsilon_{i,\Delta}^m(x)&=\left(\frac{h}{2}\right)^{i}{\cal J}_i^m(2(x-M)/h), \, i=2m+2, 2m+3, \ldots, p, \tag{38} \\ \Upsilon_{p+1,\Delta}^m(x)&=\left(\frac{h}{2}\right)^{p+1}{\cal J}_{p+1}^m(2(x-M)/h)\equiv\sum_{i=0}^{p+1}{\hat{a}}_{p+1-2i}^{p+1,m}(x-M)^{p+1-2i}, \, {\hat{a}}_{p+1}^{p+1,m}=1, \tag{39} \\ \Upsilon_{p-2m-1,\Delta}^{-m-2}(x)&=\left(\frac{h}{2}\right)^{p-2m-1}{\cal J}_{p-2m-1}^{-m-2}(2(x-M)/h), \tag{40} \end{align}\] where from 7 \[\Upsilon_{p+1,\Delta}^m(x)=(x-M+h/2)^{m+1}(x-M-h/2)^{m+1}\Upsilon_{p-2m-1,\Delta}^{-m-2}(x), \label{eq:upsilonp1a}\tag{41}\] and the \(\Upsilon_{i,\Delta}^m(x)\), \(i=2m+2, 2m+3, \ldots, p+1\) are monic. Further, let \[\begin{align} W_{p,\Delta}^m(x)&=\begin{cases} \sum_{i=0}^m(W_{i,L}\Phi_{i,L,\Delta}^m(x)+W_{i,R}\Phi_{i,R,\Delta}^m(x)) & p=2m+1, \\ \sum_{i=0}^m(W_{i,L}\Phi_{i,L,\Delta}^m(x)+W_{i,R}\Phi_{i,R,\Delta}^m(x))+\sum_{i=2m+2}^p W_{i,M}\Upsilon_{i,\Delta}^m(x) & p \ge 2m+2, \end{cases} \label{eq:interpolant} \end{align}\tag{42}\] be an interpolant on \(\Delta\) such that \[\begin{align} W_{i,L}&=u^{(i)}(M-h/2), \, i=0, 1, \ldots m, \tag{43} \\ W_{i,R}&=u^{(i)}(M+h/2), \, i=0, 1, \ldots m, \, \tag{44} \end{align}\] and, if \(p \ge 2m+2\), the \(W_{i,M}\) are chosen so that \[W_{p,\Delta}^m(M_i)=u(M_i), \, i=2, 3, \ldots, p-2m, \label{eq:middle}\tag{45}\] with \(M_i=M+h\mu_i/2\) where \(\mu_i\) are the \(p-2m-1\) roots of \(\Upsilon_{p-2m-1,\Delta}^{-m-2}(\mu)\), \(i=2, 3, \ldots, p-2m\) and \(M_1=M-h/2\) and \(M_{p-2m+1}=M+h/2\). Additionally, let \(\pi^{\Delta}\) be the linear operator that projects functions in \(C^{\infty}({\bar \Delta})\) onto the space of polynomials of degree \(p\) on \({\bar \Delta}\) according to 42 45 and let \[\begin{align} v(x)=&\sum_{i=0}^{p+s+1} \frac{(x-M)^i}{i!}u^{(i)}(M), \tag{46} \\ {\hat{L}}_{i,\Delta}^{m,p}(x)=&\frac{(x-M_1)^{m+1}}{(M_i-M_1)^{m+1}}\frac{(x-M_{p-2m+1})^{m+1}}{M_i-M_{p-2m+1})^{m+1}}\prod_{j=2,j \neq i}^{p-2m}\frac{(x-M_j)}{(M_i-M_j)}, i=2, 3, \ldots, p-2m, \tag{47} \\ {\hat{A}}_{l,\Delta}^{m,p}(x)=&\frac{\sum_{j=0}^{m-l}A_{j,l}(x-M_1)^{j+l}}{l!}\frac{(x-M_{p-2m+1})^{m+1}\prod_{j=2}^{p-2m}(x-M_j)}{(M_1-M_{p-2m+1})^{m+1}\prod_{j=2}^{p-2m}(M_1-M_j)}, \tag{48} \\ {\hat{B}}_{l,\Delta}^{m,p}(x)=&\frac{\sum_{j=0}^{m-l}B_{j,l}(x-M_{p-2m+1})^{j+l}}{l!}\frac{(x-M_1)^{m+1}\prod_{j=2}^{p-2m}(x-M_j)}{(M_{p-2m+1}-M_1)^{m+1}\prod_{j=2}^{p-2m}(M_{p-2m+1}-M_j)}, \tag{49} \\ l&=0, 1, \ldots, m, \end{align}\] where the \(A_{j,l}\) and \(B_{j,l}\) are determined by requiring that \[\begin{align} ({\hat{A}}_l^{m,p}(M_1))^{(k)} &= \begin{cases} 1 & l=k, \\ 0 & l \neq k, \, k=0, 1, \ldots, m, \end{cases} \tag{50} \\ ({\hat{B}}_l^{m,p}(M_{p-2m+1}))^{(k)}&= \begin{cases} 1 & l=k, \\ 0 & l \neq k, \, k=0, 1, \ldots, m. \end{cases} \tag{51} \end{align}\] Remark: The number of conditions 50 51 is larger than the number of coefficients \(A_{j,l}\) and \(B_{j,l}\), respectively. That unique values of the coefficients can be found while simultaneously meeting the conditions is demonstrated below.

Remark: If \(m=0\) the \({\hat{L}}_{i,\Delta}^{m,p}(x)\), \({\hat{A}}_{l,\Delta}^{m,p}(x)\), and \({\hat{B}}_{l,\Delta}^{m,p}(x)\) are the Lagrange interpolants.

Before proceeding to the proof of the theorem, four lemmas are proved that establish that the \({\hat{L}}_{i,\Delta}^{p,m}(x)\), \({\hat{A}}_{l,\Delta}^{p,m}(x)\), \({\hat{B}}_{l,\Delta}^{p,m}(x)\) form a basis for \(S^{p,\Delta}\). Thus, they constitute generalizations of the Lagrange interpolating polynomials. One further lemma demonstrates a key property of \(\Upsilon_{p+1,\Delta}^m(x)\).

To prove the comparable properties of the \({\hat{A}}_{l,\Delta}^{m,p}(x)\) and \({\hat{B}}_{l,\Delta}^{m,p}(x)\) it is helpful to decompose them into smaller components. To that end let \[{\hat{A}}_{l,\Delta}^{m,p}(x)=F_{l,\Delta}^m(x)G_{p,\Delta}^m(x), \label{eq:hataproduct}\tag{52}\] where \[\begin{align} F_{l,\Delta}^m(x)&=\sum_{j=0}^{m-l}A_{j,l}(x-M_1)^{j+l}, \tag{53} \\ G_{p,\Delta}^m(x)&={\bar G}_{p,\Delta}^m(x-M_{p-2m+1})^{m+1}\prod_{j=2}^{p-2m}(x-M_j), \tag{54} \\ {\bar G}_{p,\Delta}^m&=\frac{1}{(M_1-M_{p-2m+1})^{m+1}\prod_{j=2}^{p-2m}(M_1-M_j)}, \tag{55} \\ {\hat{G}}_{p,\Delta}^m(x)&=\frac{m+1}{x-M_{p-2m+1}}+\sum_{j=2}^{p-2m}\frac{1}{x-M_j}, \tag{56} \end{align}\] with analogous functions for \({\hat{B}}_{l,\Delta}^{m,p}.\)

Proof (of Theorem 6). Proofs for the cases \(m=0\) and \(m=1\) can be found in [24] and [25], respectively.

Consider the expansion \[\label{eq:uexpansion} u(x)-W_{p,\Delta}^m(x)=u(x)-v(x)+v(x)-\pi^{\Delta}v(x)+\pi^{\Delta}v(x)-W_{p,\Delta}^m(x)\tag{67}\] and address each pair in turn. For the first pair from 46 and Taylor’s Theorem it follows that \(\forall x \in (M-h/2,M+h/2) \, \exists \, {\hat{x}}\) s.t. \[u(x)-v(x)=\frac{u^{(p+s+2)}(M)}{(p+s+2)!}({\hat{x}}-M)^{p+s+2}.\] Moving to the second pair, using the linearity of \(\pi^{\Delta}\) and 46 \[v(x)-\pi^{\Delta}v(x)=\sum_{i=p+1}^{p+s+1}((x-M)^i-\pi^{\Delta}(x-M)^i)\frac{u^{(i)}(M)}{i!}\equiv\sum_{i=p+1}^{p+s+1}v_i(x)\frac{u^{(i)}(M)}{i!}.\] Since the \(v_i(x)\), \(i=p+1, p+2, \ldots, p+s+1\), satisfy ?? ?? it follows from Lemma 11 that \[v_i(x)=\Upsilon_{p+1,\Delta}^m(x)\Psi_{i-p-1}^m(x), \, i=p+1, p+2, \ldots, p+s+1,\] where \(\Psi_{i-p-1}^m(x)\) is a polynomial of degree \(i-p-1\). The coefficients of \(\Psi_i^m(x)\) can be derived following the same procedure for the case \(m=0\) in Theorem 3.1 in [24]. Finally, since \(\pi^{\Delta}(v(x)-W_{p,\Delta}^m(x)) \in S^{p,\Delta}\) from Lemma 10 it follows that \[\begin{align} &\pi^{\Delta}(v(x)-W_{p,\Delta}^m(x))=\sum_{i=2}^{p-2m}(v(\Gamma_i)-W_{p,\Delta}(\Gamma_i)) {\hat{L}}_i(x)\\ &+\sum_{i=1}^{m+1}(v^{(i)}(\Gamma_1)-W_{p,\Delta}^{(i)}(\Gamma_1)){\hat{A}}_i^{m,p}(x)+\sum_{i=1}^{m+1}(v^{(i)}(\Gamma_{p-2m+1})-W_{p,\Delta}^{(i)}(\Gamma_{p-2m+1})){\hat{B}}_i^{m,p}(x). \end{align}\] 0◻

Remark: It is easily verified that \(\forall m, ~\Psi_0^m(x)=1, ~\Psi_1^m(x)=x-M\).

Remark: This theorem generalizes Theorem 3.1 of [24] and Theorem 3 of [25].

Let the \(L^2\) and \(H^k\) norms over \(\Delta\) be denoted by \(\Vert \cdot \Vert_{0, \Delta}\), \(\Vert \cdot \Vert_{k,\Delta}, k=1, 2, \ldots m+1\), and their corresponding inner products by \((\cdot,\cdot)_{0,\Delta}\), \((\cdot,\cdot)_{k,\Delta}, \, k=1, 2, \ldots, m+1\), respectively. In addition, let the \(H^k\) seminorms over \(\Delta\) be denoted by \(\vert \cdot \vert_{k,\Delta}, k=1, 2, \ldots m+1\). Finally let \({\tilde{C}}(\cdot)\) and \({\bar C}(\cdot)\) denote generic constants that depend on a set of parameters such as the size of the domain, \(\vert \Delta \vert\), \(p\), and the order of the derivative \(k\).

To prove ?? it is helpful to rewrite ?? in the case \(s=0\) as \[u(x)-W_{p,\Delta}^m(x)=\Upsilon_{p+1,\Delta}^m(x)\frac{u^{(p+1)}(M)}{(p+1)!}+\Pi_{p+2,\Delta}^m(x) \equiv \Pi_{p+1,\Delta}^m(x)+\Pi_{p+2,\Delta}^m(x), \label{eq:theorem1a}\tag{68}\] where \[\begin{align} \Pi_{p+2,\Delta}^m(x,{\hat{x}})&=\frac{u^{(p+2)}(M)}{(p+2)!}({\hat{x}}-M)^{p+2}+\sum_{i=2}^{p-2m}(v(M_i)-u(M_i)) {\hat{L}}_{i,\Delta}^{m,p}(x) \nonumber \\ &+\sum_{l=0}^{m}(v^{(l)}(M_1)-u^{(l)}(M_1)){\hat{A}}_{l,\Delta}^{m,p}(x)+\sum_{l=0}^{m}(v^{(l)}(M_{p-2m+1})-u^{(l)}(M_{p-2m+1})){\hat{B}}_{l,\Delta}^{m,p}(x) \nonumber \\ &\equiv \sum_{i=1}^4\Pi_{p+2,\Delta}^{m,i}(x). \label{eq:xisum} \end{align}\tag{69}\] Equations 68 69 imply that a priori estimates of the error and its derivatives require bounds on the terms \[\int_{\Delta} (\Pi_{p+1,\Delta}^m(x))^{(k)} (\Pi_{p+2,\Delta}^m(x))^{(k)}dx, \, \int_{\Delta} ((\Pi_{p+2,\Delta}^m(x))^{(k)})^2dx \label{eq:crosstermsa}\tag{70}\] along with a formula for the leading term \[\int_{\Delta} ((\Pi_{p+1,\Delta}^m(x))^{(k)})^2dx. \label{eq:leadingterm}\tag{71}\] Applying Cauchy-Schwarz to 70 yields \[\int_{\Delta} (\Pi_{p+1,\Delta}^m(x))^{(k)} (\Pi_{p+2,\Delta}^m(x))^{(k)}dx \le \Vert \Pi_{p+1,\Delta}^m(x)\Vert_k\Vert \Pi_{p+2,\Delta}^m(x)\Vert_k. \label{eq:crosstermsb}\tag{72}\] Cauchy-Schwarz also implies that \[\int_{\Delta} ((\Pi_{p+2,\Delta}^m(x))^{(k)})^2dx \le 4 \sum_{i=1}^4 \Vert \Pi_{p+2,\Delta}^{m,i}\Vert_{k,\Delta}^2. \label{eq:crosstermsc}\tag{73}\]

It is helpful to split the proof of the theorem into a series of lemmas beginning with bounds on the coefficients of \({\hat{L}}_{i,\Delta}^{p,m}(x)\), \({\hat{A}}_{i,\Delta}^{p,m}(x)\), and \({\hat{B}}_{i,\Delta}^{p,m}(x)\). These bounds depend on the separation between the distinct roots of \(\Upsilon_{p+1,\Delta}^m(x)\), \(M_i\), \(i=1, 2, \ldots, p-2m+1\), that is, the scaled roots of \(\Upsilon_{p-2m-1,\Delta}^{-m-2}(\mu)\) along with the endpoints of \(\Delta\).

The next step is to derive bounds on \(({\hat{L}}_{i,\Delta}^{p,m}(x))^{(k)}\), \({\hat{A}}_{i,\Delta}^{p,m}(x)^{(k)}\), and \({\hat{B}}_{i,\Delta}^{p,m}(x)^{(k)}\). These are then used to obtain bounds on the terms \(({\hat{A}}_l^{m,p}(x))^{(k)}\) and \(({\hat{B}}_l^{m,p}(x))^{(k)}\) and from there bounds on norms of 72 73 .

Remark: Equations ?? ?? can be used to calculate the values of \({\cal J}_p^m(\mu)\) at quadrature points.

Proof (of Theorem 12). From 68 \[\begin{align} &(u-W_{p,\Delta}^m,u-W_{p,\Delta}^m)_{k,\Delta}^2=(\Pi_{p+1,\Delta}^m+\Pi_{p+2,\Delta}^m,\Pi_{p+1,\Delta}^m+\Pi_{p+1,\Delta}^m)_{k,\Delta} \nonumber \\ &=(\Pi_{p+1,\Delta}^m,\Pi_{p+1,\Delta}^m)_{k,\Delta}+(\Pi_{p+2,\Delta}^m,\Pi_{p+2,\Delta}^m)_{k,\Delta}+2(\Pi_{p+1,\Delta}^m,\Pi_{p+2,\Delta}^m)_{k,\Delta} \label{eq:theorem2} \end{align}\tag{103}\] Bounds on the first two terms on the right-hand side of the second step of 103 come from ?? and ?? , respectively. The bound on the third term follows from ?? and ?? together with Cauchy-Schwarz. 0◻

Remark: This theorem generalizes Corollary 2.6 of [27] and Corollary 4 of [25].

From the lemmas and Theorem 6, point-wise superconvergence can now be demonstrated.

Thus, the error of the \(k^{th}\) derivative of the Hermite-Jacobi interpolant at the roots of \((\Upsilon_{p+1,\Delta}^m(x))^{(k)}\) is one order higher than the expected pointwise error and thus exhibits pointwise superconvergence. This implies the counter-intuitive result that the number of superconvergence points increases by one for every increase in the order of the derivative.

In the case \(k=m+1\) the interpolant error simplifies considerably.

Remark: This corollary generalizes Theorem 4.2 of [27] and Corollary 5 of [25].

Consider now the domain \(\Omega\) and partition it into \(N\) intervals \(\Omega_j=(x_{j-1},x_j), j=1, 2, \ldots, N\), letting \[\begin{align} \Gamma_{\Omega} & \equiv \{ a=x_0 < x_1 < \cdots < x_N=b \}. \end{align}\] Furthermore, let \(h_j=x_j-x_{j-1}\) and \(m_j=(x_{j-1}+x_j)/2\), be the interval (element) lengths and midpoints of \(\Omega_j, j=1, 2, \ldots , N\), respectively, with \(H={\rm max}_j h_j\), \(\beta H={\rm min}_jh_j\) where \(\beta\) is bounded from below independently of \(N\). Let \(S^{p,m,\Gamma_{\Omega}}\) be the space of piecewise \(C^{m+1}(\Gamma_{\Omega})\) polynomials of degree \(p \ge 2m+1\), \(j=1, 2, \ldots, N\). Consider \(W_{p,m,\Gamma_{\Omega}} \in S^{p,m,\Gamma_{\Omega}}\) \[\begin{align} W_{p,m,\Gamma_{\Omega}}&=\begin{cases} \sum_{j=1}^{N+1}\sum_{l=0}^mW_{l,j}\Phi_{l,j}^m(x) & p=2m+1, \\ \sum_{j=1}^{N+1}\sum_{l=0}^mW_{l,j}\Phi_{l,j}^m(x)+\sum_{j=1}^N\sum_{i=2m+2}^p W_{i,j}\Upsilon_{i,j}^m(x) & p \ge 2m+2, \end{cases} \label{eq:HLOMEGAa} \end{align}\tag{104}\] where \[\Phi_{l,j}^m(x)=\begin{cases} \Phi_{l,L,m_j}^m(x) & x \in \Omega_j, \\ \Phi_{l,R,m_{j-1}}^m(x) & x \in \Omega_{j-1}, \\ 0 & {\rm otherwise}, \end{cases} \label{eq:HLOMEGAb}\tag{105}\] and \[\Upsilon_{i,j}^m(x)=\begin{cases} \Upsilon_{i,\Omega_j}^m(x) & x \in \Omega_j, \\ 0 & {\rm otherwise}. \end{cases} \label{eq:HLOMEGAc}\tag{106}\] The \(W_{i,j}\) are determined by applying 43 45 to each element \(\Omega_j, j=1, 2, \ldots , N\).

Remark: This theorem generalizes Corollary 5 of [25].

4 Asymptotic Equivalence↩︎

Let \(u(x) \in C_0^{p+2}({\bar \Omega})\) be the solution of the linear equation \[\begin{align} f(x)&=\sum_{k=1}^{m+1}((-1)^kp_k(x)u^{(k)}(x))^{(k)}+p_0(x)u(x), \, x \in \Omega, \tag{107} \\ &u^{(k)}(a)=u^{(k)}(b)=0, \, k=0, 1, \ldots, m, \tag{108} \end{align}\] where \(f(x), p_k(x) \in C^p({\bar \Omega}), \, k=0, 1, \ldots, m+1\) and \(p_{m+1}(x) > 0\) with \(P_{\rm max}=\max_{k \in \{0, 1, \ldots, m+1\},x \in \Omega} |p_k(x) |\) and \(p_{\rm max}^{\prime}=\max_{x \in \Omega}| p_{m+1}^{\prime}(x) |\). Then \(u(x)\) satisfies the weak form \[\begin{align} (f,v)&=\sum_{k=0}^{m+1}\sum_{j=1}^N(p_k(x)u(x),v(x))_{k,\Omega_j} \equiv \sum_{k=0}^{m+1}a_k(u,v) \equiv a(u,v), \, \forall v \in H_0^{p}(\Omega) \label{eq:solna} \end{align}\tag{109}\] together with 108 . The finite element approximation to \(u(x)\), \(U_{p,m,\Gamma_{\Omega}}(x) \in S_0^{p,m,\Gamma_{\Omega}}\), satisfies \[\begin{align} (f,V_{p,m,\Gamma_{\Delta}})&=\sum_{k=0}^{m+1}a_k(U_{p,m,\Gamma_{\Omega}},V_{p,m,\Gamma_{\Omega}})=a(U_{p,m,\Gamma_{\Omega}},V_{p,m,\Gamma_{\Omega}}), \, \forall V_{p,m,\Gamma_{\Omega}} \in S_0^{p,m,\Gamma_{\Omega}}. \label{eq:UFEMa} \end{align}\tag{110}\]

First it will be necessary to derive an extension of the Poincaré inequality.

Thus, \(\Vert W_{p,m,\Gamma_{\Omega}}-U_{p,m,\Gamma_{\Omega}} \Vert_{m+1}\) is one order higher than \(\Vert u-U_{p,m,\Gamma_{\Omega}}\Vert_{m+1}\) and \(\Vert u-W_{p,m,\Gamma_{\Omega}}\Vert_{m+1}\). This superconvergence property implies the asymptotic equivalence of the finite element solution and the Hermite-Jacobi interpolant.

5 Computational Results↩︎

The goal of this section is to examine the impact of \(m\) on the other parameters \(p\) and \(N\). While only one example is explored, the results are suggestive. To that end consider \[0=-(-\Delta)^{m+1}u+u+f(x), \, x \in [0,1], \label{eq:test}\tag{125}\] where the boundary conditions and \(f(x)\) are chosen so that \[u(x)={\rm tanh}10(x-1/2).\] The solution profile mimics an interior boundary layer at \(x=1/2\). A finite element code written in Python using the Hermite-Jacobi basis has been developed. It uses the package mpmath [28] which allows for arbitrary precision arthimetic. While the code therefore runs considerably slower, using mpmath addresses the valid concerns about the loss of precision for high order methods in general. For example, in [29] the authors note that in the one-dimensional case the impact of round-off accumulation for fourth-order problems can be "seriously damaging" due to larger condition numbers than for second-order equations. The mpmath parameter that controls precision, \(mp.dps\) is set to 120 (having run tests for both 80 and 160).

The relationship between the three parameters was explored through three tests. The first two considered the behavior of the error in the interpolant \[E_{I,m,p,N}\equiv\Vert u-W_{p,m,N} \Vert_{m+1}, \label{eq:interperror}\tag{126}\] (where \(\Gamma_{\Omega}\) has been replaced by the number of elements \(N\)) by examining its leading term. From ?? and Stirling’s approximation it follows that \[E_{I,m,p,N} \sim \frac{1}{2(2p-2m+1)}\left(\frac{he}{4(p-m)}\right)^{p-m}\left(\int_{\Omega}(u^{(p+1)}(x))^2dx \right)^{1/2}, \label{eq:testcurve}\tag{127}\] where \(p-m > e\) for \(m \ge 2\) since \(p \ge 2m+1\). In the first experiment \(p\) is left fixed so as to focus on the relationship between \(N\) and \(m\). Thus, in 127 , \(\int_{\Omega}(u^{(p+1)}(x))^2dx\) is held constant and for each \(m \in \{2, 3, \ldots 8 \}\) \(N\) is chosen to keep the left-hand side of 127 constant. In particular \(p=17\) and the interpolation error is fixed at \(1.5969 \times 10^{-4}\). To generate the curve specified by 127 was used. The results shown in Figure 1 indicate that to achieve the same error \(N\) must increase exponentially as a function of \(m\). The finite element errors (label values) for the same \(N\), \(m\), and \(p\) match the interpolant errors except when \(m=2, 3\) when they are significantly smaller. In these two cases, \(N=6, 10\), the interpolant and finite element solution are in the pre-asymptotic range.

In the second experiment the left-hand side of 127 is fixed and for each \(m\), \(N\) is initially set to be 160 until \(p\) is selected such that \(E_{I,m,p,N} \le 1.2669 \times 10^{-10} < E_{I,m,p-1,N}\). Then \(N\) is decreased until \(E_{I,m,p,N} \approx 1.2669 \times 10^{-10}\). The number of degrees of freedom (dof) as a function of \(m\) is displayed in Figure 2. The growth in the dof is linear as is the growth of \(p\) which satisfies \(p=2m+6\). The best fit slope for degrees of freedom for this example is equal to 225 which is not insubstantial. As in experiment one the finite element errors (label values) are close to the interpolant errors.

In experiment three the finite element error, \[E_{m,p,N}=\Vert u-U_{p,m,N}\Vert_{m+1}\] is compared to \(E_{m,p,2N}\) and to \(E_{I,m,p,N}\). To that end let \[\Theta=\frac{E_{m,p,N}}{E_{I,m,p,N}}, \, \kappa=\frac{E_{m,p,N}}{E_{m,p,2N}},\] Equation 110 is solved on uniform grids with \(\beta=1\) and thus \(N=1/H\). Table ¿tbl:tab:table1? shows the results for \(m=2, 5, 8\) and two values of \(p\) for each \(m\). The smaller value of \(p\) in each case, that is \(p=5, 11, 17\), corresponds to a nodal basis only. For all three values of \(m\), \(\Theta\) is 1 to five digits. If interior basis functions are added \(\Theta\) is slightly less than 1 until \(N\) is sufficiently large. The rate at which \(\kappa\) approaches its theoretical limiting value decreases as \(m\) increases. The a posteriori error estimation strategy that guided the \(h\)-, \(p\)-, and \(hp\)-adaptivity strategies for \(m=0\) and \(m=1\) [24] was more dependable when \(\Theta\) and \(\kappa\) were in the asymptotic range. While obviously limited, these experiments suggest that as \(m\) increases more computational resources will be required to achieve the same error and that the adaptive strategies will be less robust. In particular the tests indicated that h-refinement on its own becomes a less effective strategy as \(m\) increases.

6 Conclusions↩︎

An explicit formula for Hermite splines for any interelement continuity \(C^m\), \(m \ge 0\), was presented. These splines serve as the nodal basis functions for a finite element method for solving two-point boundary equations of arbitrary even order. The splines are complemented by interior basis functions formed from ultraspherical polynomials multiplied by binomial powers of order \(m+1\) at each element end point. These polynomials were then used to construct an interpolant of the solution. Formulas for the a priori error estimates of the interpolation error were presented and extensions of the Lagrange interpolants were introduced. The leading term in the interpolation error was proportional to the \(p+1\)st derivative of the solution multiplied by the next highest degree interior basis function. It was shown that the difference between the interpolant and the finite element solution in \(H^{m+1}\) is one order higher than the a priori errors of both which in turn implies the asymptotic equivalence of the interpolant and the finite element solution. Computational results suggest that the performance gap between hp- and h-refinement grows with \(m\).

The Hermite-Jacobi basis can be extended in a straightforward manner to tensor-product domains in two and three space dimensions (e.g., rectangles, rectangular prism). They can also be used in adaptive solvers for higher-order time-dependent problems. The connection between the nodal basis functions and those in the interior suggest that this framework can be employed in various types of spectral and spectral element methods. It is possible to show that the a posteriori error estimator developed for 2nd- [24], [27] and 4th-order [25] equations can be extended to any \(m\) and \(p\) which, in turn, can serve as the key driver in a hp-adaptive finite element solver. Issues regarding precision as \(m\) increases must also be resolved and adaptivity may be part of the solution.

Acknowledgements↩︎

The author thanks Tom Hagstrom for bringing to his attention the papers by Birkhoff, Schultz, Varga and by Spitzbart.

References↩︎