Accurate Chebyshev collocation solutions
for the biharmonic eigenproblem on a rectangle

Imre Boros$^\ast $

July 18, 2017.

$^\ast $Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy, P.O. Box 68-1, 3400 Cluj-Napoca, Romania, e-mail: boros.imre@ictp.acad.ro

We are concerned with accurate Chebyshev collocation (ChC) solutions to fourth order eigenvalue problems. We consider the 1D case as well as the 2D case. In order to improve the accuracy of computation we use the precondtitioning strategy for second order differential operator introduced by Labrosse in 2009. The fourth order differential operator is factorized as a product of second order operators. In order to asses the accuracy of our method we calculate the so called drift of the first five eigenvalues. In both cases ChC method with the considered preconditioners provides accurate eigenpairs of interest.

MSC. 65M70, 34L16

Keywords. spectral methods, Chebyshev collocation, preconditioning, fourth order eigenvalue problems

1 Introduction

The idea of this article started from the work of Owen [ 8 ] , where the author deduces asymptotic estimates for the first eigenvalue of the biharmonic operator on a rectangle.

It is well known that Chebyshev collocation (ChC) method involves full populated and bad conditioned differential matrices (see for instance Boyd [ 3 ] , Gheorghiu [ 4 ] , Trefethen [ 10 ] ). In order to surpass this major difficulty we use the preconditioning technique due to Labrosse [ 6 ] and compare the obtained results with those provided by the naive methods.

We analyze the effects of preconditioning $D^4_{ChC}$ with preconditioners for $D^2_{ChC}$ based on the identity $D^4_{ChC}=D^2_{ChC}\cdot D^2_{ChC}$. We deal with the biharmonic eigenproblem solved with Chebyshev collocation and study the accuracy and the numerical stability of the method. We also consider preconditioning of the biharmonic eigenproblem and study the effects of this procedure.

In section 2 we study a fourth order eigenproblem in one dimension used by Owen [ 8 ] in order to deduce the asymptotic estimates for the first and the third eigenvalues of the biharmonic operator.

In section 3 we solve by standard Chebyshev collocation the biharmonic eigenproblem and we also consider a preconditioned version of the problem.

The paper ends with some concluding remarks and open issues.

2 Fourth order operator in $\mathbb R$

In this section we study the eigenproblem associated to the following fourth order self-adjoint operator

\begin{equation} \label{fourth_order_1D} H_{h,\alpha }(u) := \frac{d^4u}{dx^4}-2\alpha \frac{d^2u}{dx^2}, \quad \forall u \in L^2([0,h]), \end{equation}

supplied with clamped boundary conditions.

Owen analyzed in [ 8 ] this operator in order to deduce asymptotic estimates for the eigenvalues of the biharmonic operator. We solve the eigenproblem associated to (??) with standard ChC as well as with preconditioned ChC.

Owen [ 8 ] gives an "exact" formula for the eigenvalues of $H_{1,\alpha }$. Namely, for $\alpha {\gt}0$, let $\beta {\gt}\alpha $ be the $n$-th solution of the transcendental equation

\begin{equation} \label{transcend} \cosh (\sqrt{\beta +\alpha })\cos (\sqrt{\beta -\alpha })-\tfrac {\alpha \sinh (\sqrt{\beta +\alpha })\sin (\sqrt{\beta -\alpha })}{\sqrt{\beta ^2-\alpha ^2}}=1. \end{equation}

Then the $n$-th order eigenvalue is

\begin{equation} \label{nth_eig} \mu _n(\alpha ,1)=\beta ^2-\alpha ^2. \end{equation}

For our numerical experiments we choose $\alpha =1 $. In this particular case we can obtain the first eigenvalue of the operator (??) by solving (2) with Newton’s method and (??). Thus, we get the approximate eigenvalue

\begin{equation} \label{n1sth_eig} \mu _1(1,1)=525.1393094840871. \end{equation}

We will use the above value in order to check the accuracy of the first eigenvalue which we have obtained from ChC. Actually the eigenvalue problem associated to the fourth order operator (??) has the following form

\begin{equation} \label{fourth_order_evp} \begin{cases} \tfrac {d^4u}{dx^4}-2\alpha \tfrac {d^2u}{dx^2}=\mu u, \qquad x\in [0,h]\\[0.2cm] u(0)=u’(0)=u(h)=u’(h)=0. \end{cases} \end{equation}

In order to use ChC we have to transform the problem from $[0,h]$ to $[-1,1]$, for this we use the transformation $t=\frac{2}{h}x-1, dx=\frac{h}{2}dt$. After the transformation we use the differentiation matrices from [ 11 ] and the method proposed by Hoepffner in [ 5 ] to impose the boundary conditions. In this way we get the following algebraic form

\begin{equation} \label{algebraic_evp} Au=\mu u, \end{equation}

where

\begin{equation} A = \widetilde{D}_{ChC}^4-2\alpha \widetilde{D}_{ChC}^2, \end{equation}

where in $\widetilde{D}$ the boundary conditions have been introduced. In order to solve the eigenvalue problem (??) we use the eig function from Matlab.

\includegraphics[scale=.44]{first_four-eps-converted-to-1.png} — Figure 1 First four eigenvalues of $H(1,\alpha )$
obtained by ChC for N=64.

\includegraphics[scale=.44]{eigenfunctions-eps-converted-to-1.png} — Figure 2 The first three eigenfunctions of $H(1,1)$ obtained by ChC for $N=64$.

Eigenvalues of $H(1,\alpha )$ and eigenvectors of $H(1,1)$.

First we solve the eigenvalue problem for different values of $\alpha $ and order the eigenvalues in increasing values. The first four of these eigenvalues are plotted in Figure 1. In [ 8 ] this figure is obtained by solving the implicit formula (2) and using the formula for the $n$-th eigenvalue (??).

Preconditioning the eigenproblem. The idea of preconditoners for spectral methods was introduced by Orszag (1980) [ 9 ] . Given a cut-off frequency $N$, it was proposed that the desired spectral solution $u_N(x)$ of the linear problem

\begin{equation} \mathcal{L}_{sp}[u(x)] = f(x) \end{equation}

to be obtained from a finite precision operator, $\mathcal{L}_{app}$, instead of $\mathcal{L}_{sp}$. This is carried out iteratively, for example with the Richardson iterative method.

Labrosse studies in [ 6 ] the best known lowest-order preconditioner for the $\frac{d^2}{dx^2}$ Chebyshev spectral operator, and the analysis is carried out taking as a guide the elliptic eigenvalue problem

\begin{equation} \label{eig_value_prob} \tfrac {d^2u}{dx^2} = \xi u(x),\quad \textnormal{for} \quad x \in (-1,1) \quad \textnormal{and} \quad \xi <0. \end{equation}

For any numerical approximation $u_N(x)$ of $u(x)$ we can write the spectral and local approximation of (??) in the following form

\begin{eqnarray*} \mathcal{L}_{sp}[u_N(x)] & =& \xi _{sp}u_N(x), \\ \mathcal{L}_{app}[u_N(x)] & =& \xi _{app}u_N(x), \quad x\in (-1,1). \end{eqnarray*}

The efficiency of this method (or, for that matter, of any other iterative method) depends on the choice of $\mathcal{L}_{app}$: it should lead to an easily invertible matrix and be as close as possible to $\mathcal{L}_{sp}$. The closer the condition number of the matrix representing $\mathcal{L}_{app}^{-1}\mathcal{L}_{sp}$ is to 1, the more efficient is the method.

Another concept which we study is the effect of the preconditioning on the non-normality of the differentiation matrices. For this we use a scalar measure, for a square complex matrix $A$ the non-normality ratio or the so called Henrici’s number $H (A)$ is defined as

\begin{equation} \label{Henrici} H(A):= \frac{(\epsilon (A^*A-AA^*))^{1/2}}{\epsilon (A)} \end{equation}

where $A^*$ is the conjugate transpose of $A$ and $\epsilon (A)$ stands for the Frobenius norm of $A$. Gheorghiu [ 4 ] deduced the estimation

\[ 0\leq H(A) \leq 2^{1/4}\approx 1.1892, \]

where $H(A)=0$ if and only if $A$ is normal.

In Table 1 and Table 2 can be observed the effect of preconditioning on the condition number and on the non-normality ratio for the second and fourth order differentiation matrices for $N =64$ and $N=100$.

	$D^2_{ChC}$	$D^{2,precond}_{ChC}$	$D^4_{ChC}$	$D^{4,precond}_{ChC}$
Henrici	0.3013	0.1469	0.3227	0.9320
Cond. no.	$7.2254 \cdot 10^4$	3.0438	$9.2812\cdot 10^8$	$1.1353\cdot 10^5$

Table 1 The Henrici and the condition numbers for second and fourth order differentiation matrices when $N = 64$.

	$D^2_{ChC}$	$D^{2,precond}_{ChC}$	$D^4_{ChC}$	$D^{4,precond}_{ChC}$
Henrici	0.3023	0.1423	0.3224	0.9647
Cond. no.	$4.4011 \cdot 10^5$	3.5175	$3.4342 \cdot 10^{10}$	$8.6004 \cdot 10^5$

Table 2 The Henrici and the condition numbers for second and fourth order differentiation matrices when $N=100$.

From the above tables we can observe that by preconditioning we get better condition numbers for the second order case and also for the fourth order differentiation matrix. With respect to the non-normality ratio it is improved only for the second order differentiation matrices. For the fourth order case the matrices get more non-normal after preconditioning. The explanation of this aspect remains an open issue.

Accuracy of the method. Now we study the accuracy of the ChC method to obtain the first eigenvalue and compare our results to the preconditioned version of our eigenvalue problem.

In Table 3 we compare the values of the first eigenvalue $\mu _1$ obtained from (??), i.e., the value (??) with $\mu _1^N$ obtained with ChC standard and $\mu _{1,precond}^N$ provided by preconditioned ChC. An important result of this paper is that by preconditioning we improve the accuracy by at least one decimal digit.

N	$\|\mu _1 - \mu _1^{N} \|$	$\|\mu _1 - \mu _{1,precond}^{N} \|$
60	$6.089 \cdot 10^{-7}$	$4.326 \cdot 10^{-9}$
80	$1.468 \cdot 10^{-6}$	$1.139 \cdot 10^{-8}$
100	$2.097 \cdot 10^{-7}$	$1.668 \cdot 10^{-8}$
120	$1.593 \cdot 10^{-7}$	$6.256 \cdot 10^{-8}$
140	$4.117 \cdot 10^{-7}$	$6.011 \cdot 10^{-8}$
160	$1.394 \cdot 10^{-6}$	$1.507 \cdot 10^{-8}$
180	$4.162 \cdot 10^{-6}$	$5.251 \cdot 10^{-9}$
200	$1.922 \cdot 10^{-6}$	$4.972 \cdot 10^{-7}$

Table 3 Errors in the first eigenvalue for $N = 60, 80,...,200$.

Numerical stability. In order to asses the accuracy and the numerical stability of the algorithm, we use the concept of the eigenvalue drift introduced by Boyd in [ 3 , p. 138 ] . The eigenvalue drift is defined as the differences between the same eigenvalues which were obtained for different values of $N$, i.e. $\delta _j:=|\lambda _j^{N_1}-\lambda _j^{N_2}|$ is the drift of the $j$-th eigenvalue computed with $N_1$ respectively $N_2$ orders of approximation. In Figure 3 we plot the drift for the standard ChC and also for the preconditioned problem, for $N=100, 120$ resp. $N = 100, 140$.

\includegraphics[scale=0.87]{drift_1D-eps-converted-to-1.png} — Figure 3 The drift for $N=100,120$ resp. for $N=100,140$.

Alternatively, there exist some theoretical results concerning the accuracy in computing the eigenvalues of the non-normal matrices. In [ 7 , p. 474 ] the authors provide some upper bounds for the distance between the exact eigenvalues and their numerical approximation. The right hand side of this inequality depends on the departure from normality of the matrix at hand. We have found this upper bound quasi uncomputable. Instead, we have provided in Table 1 and 2 the Henrici’s number which quantity the above departure. The lower values of these numbers for preconditioned $D^2_{ChC}$ matrices partially explain the superiority numerical results of the preconditioned variant.

3 The biharmonic operator in $\mathbb R^2$

For the biharmonic operator in $\mathbb {R}^2$,

\begin{equation} \label{biharmonic_operator} \Delta ^2 = u_{xxxx}+2u_{xxyy}+u_{yyyy}, \end{equation}

we consider the "clamped" plate eigenvalue problem on a rectangle,

\begin{equation} \label{biharmonic} \begin{cases} \Delta ^2 u = \lambda u, \quad \text{in} \quad [0,h]\times [0,1],\\ u =\tfrac {\partial u}{\partial n}=0 \quad \text{on} \quad \partial [0,h]\times [0,1]. \end{cases} \end{equation}

Working on the square of side length $h=1$, Bjørstad and Tjøstheim in [ 2 ] obtained the lowest eigenvalue

\begin{equation} \label{exact_eig_biharm} \lambda _1 = 1294.9339795917128081703026479744. \end{equation}

Which we will use as an "exact" value to check the accuracy of our method.

For such a two dimensional problem we naturally set up a grid based on Chebyshev points independently in each direction, called a tensor product grid. The easiest way to solve a problem on a tensor product spectral grid is to use tensor products in linear algebra, also known as Kronecker products. The discretization of our eigenproblem (14) is

\begin{equation} A=\widetilde{D}^4_{ChC}\otimes I + 2(\widetilde{D}^2_{ChC}\otimes I) \cdot (I\otimes \widetilde{D}^2_{ChC}) + I \otimes \widetilde{D}^4_ChC, \end{equation}

where in $\widetilde{D}$ the boundary conditions have been introduced.

Preconditioning. In order to carry on the preconditioning of fourth order differential matrix we use the identity (see [ 11 ] ):

\begin{equation} \label{rel_D4_D2} D^4_{ChC} = D_{ChC}^2 \cdot D_{ChC}^2. \end{equation}

From relation (??) we get the idea that preconditioning $D^4_{ChC}$ with a preconditioner for $D^2_{ChC}$ will improve the condition number and the stability of our algorithm. In Table 1 and Table 2 it is confirmed that the condition number is better with preconditioning. But in this way we get the non-normality worse than without preconditioning.

Accuracy of the method. We use use the eigenvalue obtained by
Bjørstad and Tjøstheim in [ 2 ] to check the accuracy of the method.

N	$\|\lambda _1 - \lambda _1^{N} \|$	$\|\lambda _1 - \lambda _{1,prec}^{N} \|$
10	$0.4894 \cdot 10^{0}$	$4.893873 \cdot 10^{-1}$
16	$7.1080 \cdot 10^{-4}$	$7.1080\cdot 10^{-4}$
32	$1.6667 \cdot 10^{-7}$	$1.7370\cdot 10^{-7}$

Table 4 Errors in the first eigenvalue for $N = 10, 16, 32$.

In Table 4 we see that the preconditioning does not help to obtain better accuracy, in fact the first eigenvalue obtained from the preconditioned form gives poorer accuracy then the first eigenvalue obtained from the standard discretization.

As Weideman pointed out in [ 11 ] , a disadvantage of this approach for solving two-dimensional problems is that the order of the matrices are $N^2 \times N^2$, which quickly becomes large. This is why we used only $N=32$ in Table 4 for the two dimensional biharmonic operator.

Numerical stability. In Figure 4 we have plotted the drift for the first five eigenvalues of the biharmonic operator. We observe that preconditioning makes no difference in this case.

\includegraphics[scale=0.87]{drift_2D-eps-converted-to-1.png} — Figure 4 The drift for $N=25,30$ resp. for $N=20,30$.

\includegraphics[scale=0.8]{fcgltran_coef-eps-converted-to-1.png} — Figure 5 The coefficients of the eigenfunctions expansions $u_1,\dots , u_6$ in increasing order, for $N=32$

It is well known that ChC method works in physical space. We use the Fast Fourier Transform (FFT) to compute the coefficients of expansion as Chebyshev series of the first six eigenvectors. In other words, with this transform we move from the physical space to the coefficient space. The behavior of these coefficients, i.e., the way they decrease is a fairly important clue for the accuracy of our computations.

We included in Figure 5 a semilogy plot of the eigenfunctions expansions $u_1,\dots , u_5$ in increasing order for $N=32$.

As our main result we point out that the first eigenvector is computed with best accuracy and the expansion coefficients decrease fairly smooth to values under $10^{-7}$. The worst situation is, as we have expected, with the sixth order. However, this picture shows a reliable behavior of eigenvector computation of biharmonic problem (??)-(14).

4 Concluding remarks

For the fourth order operator we obtained high accuracy for the first eigenvalue and in the case of preconditioning we increased the accuracy at least one decimal digit and decreased significantly the condition number.

In the case of the biharmonic operator, the preconditioning with the preconditioners from [ 6 ] does not makes any difference. Our intuitive idea was that $D^4_{ChC} = D^2_{ChC} \cdot D^2_{ChC}$ and if we apply the preconditioner for $D^2_{ChC}$ this will improve the results for $D^4_{ChC}$. This is partially confirmed. In Table 1 and in Table 2 we can see how the condition number is decreased after preconditioning $D^4_{ChC}$ but this does not help to obtain better accuracy (see Table 4) or better numerical stability (see Figure 4).

In a future work we propose to use preconditioners for $D^4_{ChC}$ and compare the effects of them with the results obtained in this article.

Bibliography

1: L. Bauer, E.L. Reiss, Block five diagonal matrices and the fast numerical solution of the biharmonic equation, Math. Comp., 26 (1972) no. 118, 311–326. $\includegraphics[scale=0.1]{ext-link.png}$
2: P.E. Bjørstad, P. Tjøstheim, A note on high precision solutions of two fourth order eigenvalue problems, 1998. $\includegraphics[scale=0.1]{ext-link.png}$
3: J.P. Boyd, Chebyshev and Fourier Spectral Methods, Courier Corporation, New York, 2000. $\includegraphics[scale=0.1]{ext-link.png}$
4: C.I. Gheorghiu, Spectral Methods for Non-Standard Eigenvalue Problems: Fluid and Structural Mechanics and Beyond, Springer Science & Business, 2014. $\includegraphics[scale=0.1]{ext-link.png}$
5: J. Hoepffner, Implementation of boundary conditions, 2007, (accessed on August 2, 2017 http://www.lmm.jussieu.fr/ hoepffner/boundarycondition.pdf). $\includegraphics[scale=0.1]{ext-link.png}$
6: G. Labrosse, The piecewise-linear Finite Volume scheme: The best known lowest-order preconditioner for the $\frac{d^2}{dx^2}$ Chebyshev spectral operator, J. Comp. Phys., 228 (2009) no. 12, 4491–4509. $\includegraphics[scale=0.1]{ext-link.png}$
7: W. Kahan, B. Parlett, E. Jiang, Residual bounds on approximate eigensystems of nonnormal matrices, SIAM J. Numer. Anal., 19 (1982) no. 3, 470–484. $\includegraphics[scale=0.1]{ext-link.png}$
8: M.P. Owen, Asymptotic first eigenvalue estimates for the biharmonic operator on a rectangle, J. Differential Equations, 190 (1997) no. 1, 166–190. $\includegraphics[scale=0.1]{ext-link.png}$
9: S.A. Orszag, Spectral methods for problems in complex geometries, J. Comp. Phys, 37 (1980) no. 1, 70–92. $\includegraphics[scale=0.1]{ext-link.png}$
10: L.N. Trefethen, Spectral Methods in Matlab, SIAM, Philadelphia, PA, 2000. $\includegraphics[scale=0.1]{ext-link.png}$
11: J.A.C. Weideman, S.C. Reddy, A MATLAB differentiation matrix suite, ACM Trans. Math. Software, 26 (2000) no. 4, 465-519. $\includegraphics[scale=0.1]{ext-link.png}$

This work is licensed under a Creative Commons Attribution 4.0 International License.

N	\(\|\mu _1 - \mu _1^{N} \|\)	\(\|\mu _1 - \mu _{1,precond}^{N} \|\)
60	\(6.089 \cdot 10^{-7}\)	\(4.326 \cdot 10^{-9}\)
80	\(1.468 \cdot 10^{-6}\)	\(1.139 \cdot 10^{-8}\)
100	\(2.097 \cdot 10^{-7}\)	\(1.668 \cdot 10^{-8}\)
120	\(1.593 \cdot 10^{-7}\)	\(6.256 \cdot 10^{-8}\)
140	\(4.117 \cdot 10^{-7}\)	\(6.011 \cdot 10^{-8}\)
160	\(1.394 \cdot 10^{-6}\)	\(1.507 \cdot 10^{-8}\)
180	\(4.162 \cdot 10^{-6}\)	\(5.251 \cdot 10^{-9}\)
200	\(1.922 \cdot 10^{-6}\)	\(4.972 \cdot 10^{-7}\)

N	\(\|\lambda _1 - \lambda _1^{N} \|\)	\(\|\lambda _1 - \lambda _{1,prec}^{N} \|\)
10	\(0.4894 \cdot 10^{0}\)	\(4.893873 \cdot 10^{-1}\)
16	\(7.1080 \cdot 10^{-4}\)	\(7.1080\cdot 10^{-4}\)
32	\(1.6667 \cdot 10^{-7}\)	\(1.7370\cdot 10^{-7}\)

Accurate Chebyshev collocation solutions
for the biharmonic eigenproblem on a rectangle

1 Introduction

2 Fourth order operator in \(\mathbb R\)

3 The biharmonic operator in \(\mathbb R^2\)

4 Concluding remarks

Bibliography

	\(D^2_{ChC}\)	\(D^{2,precond}_{ChC}\)	\(D^4_{ChC}\)	\(D^{4,precond}_{ChC}\)
Henrici	0.3013	0.1469	0.3227	0.9320
Cond. no.	\(7.2254 \cdot 10^4\)	3.0438	\(9.2812\cdot 10^8\)	\(1.1353\cdot 10^5\)

Accurate Chebyshev collocation solutions for the biharmonic eigenproblem on a rectangle

1 Introduction

2 Fourth order operator in \(\mathbb R\)

3 The biharmonic operator in \(\mathbb R^2\)

4 Concluding remarks

Bibliography

Accurate Chebyshev collocation solutions
for the biharmonic eigenproblem on a rectangle