Almost Exact Computation of Eigenvalues in Approximate Differential Problems

Abstract

Differential eigenvalue problems arise in many fields of Mathematics and Physics, often arriving, as auxiliary problems, when solving partial differential equations. In this work, we present a method for eigenvalues computation following the Tau method philosophy and using Tau Toolbox tools. This Matlab toolbox was recently presented and here we explore its potential use and suitability for this problem. The first step is to translate the eigenvalue differential problem into an algebraic approximated eigenvalues problem. In a second step, making use of symbolic computations, we arrive at the exact polynomial expression of the determinant of the algebraic problem matrix, allowing us to get high accuracy approximations of differential eigenvalues.

1. Introduction

Finding eigenfunctions of differential problems can be a hard task, at least for some classical problems. Among others, we can find literature in Sturm–Liouville problems, in Mathieu problems or in Orr–Sommerfeld problems describing the difficulties involved in the resolution of those problems [1,2,3,4,5,6,7,8]. The first difficulty consists of finding accurate numerical approximations for the respective eigenvalues.

In this work, we present a procedure based on the Ortiz and Samara’s operational approach to the Tau method described in [9], where the differential problem is translated into an algebraic problem. This is achieved using the called operational matrices that represent the action of differential operators in a function. We have deduced explicit formulae for the elements of these matrices [10,11] obtained by performing operations on the bases of orthogonal polynomials and, for some families, we have exact formulae, which enables the construction of very accurate operational matrices. The Tau method has already been used for these kinds of problems [5,9,12,13]; however, our work on matrix calculation formulas adds efficiency and precision to the method.

Our main purpose is to use the Tau Toolbox, a Matlab numerical library that is being developed by our research group [14,15,16]. This library allows a stable implementation of the Tau method for the construction of accurate approximate solutions for integro-differential problems. In particular, the construction of the operational matrices is done automatically. These facts led us to think that the Tau Toolbox seems to be useful for these kinds of problems.

Finally, operating with symbolic variables, we define the determinant of those matrices as polynomials and use its roots as eigenvalues’ approximations.

We present some examples showing that, using this technique in the Tau Toolbox, we are able to obtain results comparable with those reported in the literature and sometimes even better.

2. The Tau Method

Let $\mathcal{D}:\mathbb{E}\mapsto \mathbb{F}$ be an order $\nu$ differential operator, where $\mathbb{E}$ and $\mathbb{F}$ are some function spaces, and let $g_i:\mathbb{E}\mapsto \mathbb{R},i=1,\dots ,\nu$ be $\nu$ functionals representing boundary conditions, so that

$$\left\{\begin{array}{c}\mathcal{D}y=f,f\in \mathbb{F},\hfill \\ g_i(y)={\varphi }_i,i=1,\dots ,\nu \hfill \end{array}\right.$$

(1)

is a well posed differential problem.

2.1. The Tau Method Principle

A particular implementation of the Tau method depends on the choice of an orthogonal basis for $\mathbb{F}$. A sequence of orthogonal polynomials ${\{P_n(x)\}}_{n=0}^{\infty }$ with respect to the weight function $w(x)$ on a given interval of orthogonality $[a,b]$ satisfies

$$\langle P_i,P_j\rangle ={\int }_a^{b}w(x)P_i(x)P_j(x)dx=w_i{\delta }_{ij},$$

where $w_i=\langle P_i,P_i\rangle$ and ${\delta }_{ij}$ is the Kronecker delta [17].

Let $\mathbb{P}$ be the space of algebraic polynomials of any degree and let us suppose that $\mathbb{P}$ is dense in $\mathbb{F}$; then, the solution y of (1) has a series representation $y\sim {\sum }_{j\ge 1}{a}_j{P}_{j-1}$. A polynomial approximation of degree $n-1\in \mathbb{N}$ is achieved by

$$y_n=\sum _{j=1}^n{a}_{n,j}{P}_{j-1}={\mathcal{P}}_n{\mathsf{a}}_n,$$

(2)

where ${\mathcal{P}}_n=[P_0,P_1,\dots ,P_{n-1}]$ and ${\mathsf{a}}_n={[a_{n,1},\dots ,{\mathsf{a}}_{n,n}]}^T$.

In the Tau method sense, $y_n$ is a polynomial satisfying the boundary conditions in (1) and solving the differential equation with a residual ${\tau }_n=\mathcal{D}{y}_n-f$ of maximal order. Thus, the differential problem is reduced to an algebraic one of finding the n coefficients $a_{n,j},j=1,\dots ,n$ in (2) such that

$$\left\{\begin{array}{c}{g}_i(y_n)={\varphi }_i,i=1,\dots ,\nu ,\hfill \\ \langle P_{i-1},\mathcal{D}{y}_n-f\rangle =0,i=1,\dots ,n-\nu ,\hfill \end{array}\right.$$

(3)

and so the residual ${\tau }_n=\mathcal{O}(P_{n-\nu })$.

2.2. Operational Formulation

For a given $n\in \mathbb{N},n>\nu$, we define the matrix

$${\mathsf{T}}_n=\left[\begin{array}{c}{\mathsf{B}}_{\nu \times n}\hfill \\ {\mathsf{D}}_{(n-\nu )\times n}\hfill \end{array}\right]={(t_{i,j})}_{i,j=1}^n:t_{i,j}=\left\{\begin{array}{c}{g}_i(P_{j-1}),i=1,\dots ,\nu \hfill \\ {\displaystyle \frac{\langle P_{i-\nu -1},\mathcal{D}{P}_{j-1}\rangle }{w_{i-\nu -1}},i=\nu +1,\dots ,n}\hfill \end{array}\right.$$

(4)

and the vector

$${\mathsf{b}}_n=\left[\begin{array}{c}{\mathsf{g}}_{\nu }\hfill \\ {\mathsf{f}}_{n-\nu }\hfill \end{array}\right]={(b_i)}_{i=1}^n:b_i=\left\{\begin{array}{c}{\varphi }_i,i=1,\dots ,\nu \hfill \\ {\displaystyle \frac{\langle P_{i-\nu -1},f\rangle }{w_{i-\nu -1}},i=\nu +1,\dots ,n.}\hfill \end{array}\right.$$

If, in problem (1), the differential operator $\mathcal{D}$ is linear and the $g_j$ are $\nu$ linear functionals, then problem (3) can be put in matrix form as

$${\mathsf{T}}_n{\mathsf{a}}_n={\mathsf{b}}_n.$$

The matrix ${\mathsf{T}}_n$, called the Tau matrix, can be evaluated from operational matrices, that is, matrices translating into coefficients vectors the action of a differential operator $\mathcal{D}$ in a function y.

Proposition 1.

Let $\mathcal{P}=[P_0,P_1,\dots ]$ be an orthogonal polynomial basis, $y=\mathcal{P}\mathsf{a}$ and $\mathsf{M},\mathsf{N}$ infinite matrices such that

$$x\mathcal{P}=\mathcal{P}\mathsf{M}\mathit{and}{\displaystyle \frac{d}{dx}}\mathcal{P}=\mathcal{P}\mathsf{N}.$$

Then, for each $k\in {\mathbb{N}}_0$,

$$x^{k}y=\mathcal{P}{\mathsf{M}}^k\mathsf{a}\mathit{and}{\displaystyle \frac{d^k}{d{x}^k}}y=\mathcal{P}{\mathsf{N}}^k\mathsf{a}.$$

(5)

Proof. 

For $k=1$, the result is true by hypothesis. Now, supposing that (5) is true for a $k\in \mathbb{N}$, then

$$x^{k+1}y=x(x^{k}y)=(x\mathcal{P}){\mathsf{M}}^k\mathsf{a}=\mathcal{P}{\mathsf{M}}^{k+1}\mathsf{a}$$

and

$${\displaystyle \frac{d^{k+1}}{d{x}^{k+1}}}y={\displaystyle \frac{d}{dx}}({\displaystyle \frac{d^k}{d{x}^k}}y)=({\displaystyle \frac{d}{dx}}\mathcal{P}){\mathsf{N}}^k\mathsf{a}=\mathcal{P}{\mathsf{N}}^{k+1}\mathsf{a},$$

ending the proof by induction. □

The following result generalizes the algebraic representation from the previous proposition to differential operators.

Corollary 1.

Let $\mathcal{D}:{\mathbb{P}}_n\mapsto {\mathbb{P}}_{n+h}$ be a linear differential operator with polynomial coefficients

$$\mathcal{D}=\sum _{i=0}^{\nu }{p}_i{\displaystyle \frac{d^i}{d{x}^i}},p_i\in {\mathbb{P}}_{n_i}$$

(6)

and let $h={max}_{i=0,\dots ,\nu }\{n_i-i\}$.

If $y_n={\mathcal{P}}_n{\mathsf{a}}_n$, then $\mathcal{D}{y}_n={\mathcal{P}}_{n+h}{\mathsf{D}}_{(n+h)\times n}{\mathsf{a}}_n$ with

$${\mathsf{D}}_{(n+h)\times n}=\sum _{i=0}^{\nu }{p}_i(\mathsf{M}){\mathsf{N}}_{(n-i)\times n}^i,$$

where ${\displaystyle p_i(\mathsf{M})=\sum _{k=0}^{n_i}{p}_{ik}{\mathsf{M}}_{(n+h)\times (n-i)}^k}$ when ${\displaystyle p_i=\sum _{k=0}^{n_i}{p}_{ik}{x}^k}$, and ${\mathsf{A}}_{m\times n}^k$ denotes the main $m\times n$ block of the matrix ${({\mathsf{A}}_p)}^k$ with $p=max\{m,n\}$.

In [9], the authors discussed the application of this operational formulation of the Tau method to the numerical approximation of eigenvalues defined by differential equations. They proved that, for a differential eigenvalue problem, where in (1)

$$D=\sum _{r=0}^t{\lambda }^r{D}_r$$

and $\lambda$ is a parameter, the zeros of $det({\mathsf{T}}_n(\lambda ))$ approach the eigenvalues of (1).

2.3. Tau Matrices’ Properties

Given that we are dealing with a general orthogonal polynomial basis, instead of particular cases like Chebyshev or Legendre, we can only make assumptions about general properties of Tau matrices $T_n$. Anyway, we can’t expect to have symmetric matrices and, in general, they can be considered sparse but with a low level of sparsity.

Since $\mathcal{P}$ in Proposition 1 is an orthogonal basis, then $\mathsf{M}$ is the tridiagonal matrix with the coefficients of its three term recurrence relation. Therefore, for problems with polynomial coefficients, matrices $p_i(\mathsf{M})$ of Corollary 1 are banded matrices, with all non-zero elements between the $\pm n_i$ diagonals.

Matrices $\mathsf{N}$ are always strictly upper triangular and so $p_i(M)N^i$ are $n_i-i$ upper Hessenberg matrices. The resulting ${\mathcal{D}}_{(n-\nu )\times n}$ block of $T_n$ defined in (4) is a general h upper Hessenberg matrix.

Moreover, one advantage of the Tau method is its ability to deal with boundary conditions, allowing the treatment of any linear combination of values of y and of its derivatives for $g_i$ in (1). Thus, the $\nu \times n$ block $B_{\nu \times n}$ in $T_n$ is usually dense, with its entries $g_i(P_j)$, made by linear combinations of orthogonal polynomial values $P_j(x_k)$ and of its derivatives $P_j^{(l)}(x_k)$, in prescribed abscissas $x_k$.

Assembling those blocks $B_{\nu \times n}$ and ${\mathcal{D}}_{(n-\nu )\times n}$ in $T_n$, we get an $\nu +h$ upper Hessenberg matrix.

For some problems, whose dependence on the eigenvalues $\lambda$ is verified only in the differential equation, we can use Schur complements to reduce matrix sizes. Considering matrix $T_n$ in (4) partitioned as

$${\mathsf{T}}_n=\left[\begin{array}{cc}{\mathsf{B}}_1\hfill & {\mathsf{B}}_2\hfill \\ {\mathsf{D}}_1\hfill & {\mathsf{D}}_2,\hfill \end{array}\right]$$

where ${\mathsf{B}}_1$ is $\nu \times \nu$ and the other blocks are partitioned accordingly. If ${\mathsf{B}}_1$ is non-singular, then

$$det({\mathsf{T}}_n)=det({\mathsf{B}}_1)det({\mathsf{D}}_2-{\mathsf{D}}_1{\mathsf{B}}_1^{-1}{\mathsf{B}}_2)$$

and the problem is reduced to solve

$$det({\mathsf{C}}_n)=0,{\mathsf{C}}_n={\mathsf{D}}_2-{\mathsf{D}}_1{\mathsf{B}}_1^{-1}{\mathsf{B}}_2,$$

(7)

reducing to $n-\nu$ the problem dimension. In the worst case, when $B_1$ is singular, we have to work with the $n\times n$ matrix $T_n$.

In the following sections, we illustrate the application of the Tau method to approximate eigenvalues in some classical problems.

3. Problems with Polynomial Coefficients

Sturm–Liouville problems arise from vibration problems in continuum mechanics. The general form of a fourth order Sturm–Liouville equation is

$${(p(x)y^{″}(x))}^{″}-{(q(x)y^{\prime }(x))}^{\prime }+r(x)y(x)=\lambda \mu (x)y(x),a<x<b$$

(8)

with appropriate initial and boundary conditions, where $a<b\in \mathbb{R}$, $p,q,r$, and $\mu$ are given piecewise continuous functions, with $p(x)>0$ and $\mu (x)>0$. These conditions mean that (8) has an infinite sequence of real eigenvalues, bounded from above, and each one has multiplicity of at most 2 [1].

If p and q are differentiable functions, it is an elementary task to give (8) the form

$$p(x)y^{(4)}(x)+2p^{\prime }(x)y^{‴}(x)+(p^{″}(x)-q(x))y^{″}(x)-q^{\prime }(x)y^{\prime }(x)+(r(x)-\lambda \mu (x))y(x)=0.$$

From this equation, we derive the operational matrix for the general form of the fourth order Sturm–Liouville differential operator associated with (8)

$$\mathsf{D}=p(\mathsf{M}){\mathsf{N}}^4+2p^{\prime }(\mathsf{M}){\mathsf{N}}^3+(p^{″}(\mathsf{M})-q(\mathsf{M})){\mathsf{N}}^2-q^{\prime }(\mathsf{M})\mathsf{N}+r(\mathsf{M})-\lambda \mu (\mathsf{M}).$$

(9)

Assuming that coefficients $p,q,r$, and $\mu$ are polynomials, or convenient polynomial approximations of the coefficient functions, then the height of this differential operator is well defined as

$$h=max\{deg(p)-4,deg(q)-2,deg(r),deg(\mu )\},$$

where $deg(.)$ is the polynomial degree. One consequence of Corollary 1 is that to evaluate the block ${\mathsf{D}}_{(n-\nu )\times n}$ in (4) we have to apply (9) with $\mathsf{M}$ and $\mathsf{N}$ truncated to its first $n+h$ lines and columns.

The Tau matrix of a fourth order Sturm–Liouville problem is the $n\times n$ matrix $T_n=[{\mathsf{B}}_{4\times n};{\mathsf{D}}_{(n-4)\times n}]$, where ${\mathsf{B}}_{4\times n}$ is the $4\times n$ matrix representing boundary conditions and ${\mathsf{D}}_{(n-4)\times n}$ is the first $(n-4)\times n$ main block of $\mathsf{D}$.

Example 1.

Consider the Sturm–Liouville boundary value problem

$$\left\{\begin{array}{c}{y}^{(4)}(x)=\lambda y(x),0<x<1,\hfill \\ y(0)=y(1)=y^{\prime }(0)=y^{″}(1)=0,\hfill \end{array}\right.$$

(10)

whose exact eigenvalues satisfy [1,2]

$$tanh(\sqrt[4]{\lambda })-tan(\sqrt[4]{\lambda })=0.$$

(11)

In that case $\mathsf{D}={\mathsf{N}}^4-\lambda I$, where I is the identity matrix, and the boundary conditions can be represented by ${\mathsf{B}}_n=[v_0;v_1;v_0\mathsf{N};v_1{\mathsf{N}}^2]$, where $v_0=[1,-1,\cdots ,{(-1)}^{n-1}]$ and $v_1=[1,1,\cdots ,1]$ are length n line vectors with the polynomial base values in the boundary domain.

For each $n>4$, ${\mathsf{C}}_n$ in (7) is an $n-4$ square matrix and its determinant an $n-4$ degree polynomial. We use the Matlab function roots to find its zeros and we inspect their accuracy by testing if they satisfy relation (11).

In Figure 1, we present $|tanh(\sqrt[4]{{\lambda }_{n,k}})-tan(\sqrt[4]{{\lambda }_{n,k}})|$ for $k=1,\dots ,10$ the first 10 eigenvalues approximations obtained with $n=21,28,35$ and $n=42$, with Chebyshev and Legendre bases shifted to $[0,1]$.

Example 2.

A very similar problem, presented as the clamped rod problem in [12,18], is

$$\left\{\begin{array}{c}{y}^{(4)}(x)=\lambda y(x),-1<x<1,\hfill \\ y(\pm 1)=y^{\prime }(\pm 1)=0.\hfill \end{array}\right.$$

(12)

In that case, and whenever we have a symmetric problem and a symmetric base, the matrix $C_n$ in (7) has zeros intercalating all non-zero elements. We can reduce the problem dimension, defining two matrices $C_O$ and $C_E$ with respectively the odd and the even entries of $C_n$, then $det(C_n)=det(C_O)det(C_E)$. In that case, since $C_n$ is an 4-upper Hessenberg matrix, $C_O$ and $C_E$ are 2-upper Hessenberg matrices. The sparsity pattern of those two matrices, in Legendre basis and with $n=52$, are showed in Figure 2.

The first 14 eigenvalues, evaluated with an $16\times 16$ matrix, are presented in [12]. In

Table 1. Eigenvalues of Example 2 presented in [12] and ${\lambda }_{n,k}$ with $n=20$ and $n=52$ in Legendre basis.

k	${\mathit{\lambda }}_{\mathit{k}}$ [12]	${\mathit{\lambda }}_{16,\mathit{k}}$	${\mathit{\lambda }}_{52,\mathit{k}}$	${\mathit{\lambda }}_{52,\mathit{k}}/{\mathit{k}}^4$
1	3.12852439 $\times {10}^1$	3.128524385877707 $\times {10}^1$	3.128524385877221 $\times {10}^1$	31
2	2.377210675 $\times {10}^2$	2.377210675311160 $\times {10}^2$	2.377210675311198 $\times {10}^2$	15
3	9.136018866 $\times {10}^2$	9.136018831954221 $\times {10}^2$	9.136018831951466 $\times {10}^2$	11
4	2.4964874758 $\times {10}^3$	2.496487437857343 $\times {10}^3$	2.496487437856835 $\times {10}^3$	9.8
5	5.5710074688 $\times {10}^3$	5.570962978086419 $\times {10}^3$	5.570962978573987 $\times {10}^3$	8.9
6	1.08631975968 $\times {10}^4$	1.086758221396396 $\times {10}^4$	1.086758221697812 $\times {10}^4$	8.4
7	1.93928004466 $\times {10}^4$	1.926303581823010 $\times {10}^4$	1.926302825661405 $\times {10}^4$	8.0
8	3.05369477203 $\times {10}^4$	3.178016789424974 $\times {10}^4$	3.178009645408997 $\times {10}^4$	7.8
9	6.03075735206 $\times {10}^4$	4.960468438481630 $\times {10}^4$	4.958769590877672 $\times {10}^4$	7.6
10	7.11035649235 $\times {10}^4$	7.407231618213559 $\times {10}^4$	7.400084934912209 $\times {10}^4$	7.4
11	3.597677558196 $\times {10}^5$	1.091110089932048 $\times {10}^5$	1.064806931408837 $\times {10}^5$	7.3
12	3.856105241227 $\times {10}^5$	1.548144789089380 $\times {10}^5$	1.486344772858071 $\times {10}^5$	7.2
13	1.62401642422808 $\times {10}^7$	2.619042084314734 $\times {10}^5$	2.022155654215451 $\times {10}^5$	7.1
14	1.71337968904269 $\times {10}^7$	3.728032882888538 $\times {10}^5$	2.691234348268295 $\times {10}^5$	7.0

, we compare those values with our results in the Legendre basis and with $n=16$ and $n=52$. We present values of ${\lambda }_{52,k}/k^4$, which allows us to verify that our estimates satisfy the property that the kth eigenvalue is proportional to $k^4$ [18].

Example 3.

Consider the following Sturm–Liouville problem with non-null q and r coefficients

$$\left\{\begin{array}{c}{y}^{(4)}(x)-{(\alpha x^2{y}^{\prime }(x))}^{\prime }+(\beta x^4-\alpha )y(x)=\lambda y(x),0<x<5,\hfill \\ y(0)=y^{″}(0)=y(5)=y^{″}(5)=0,\hfill \end{array}\right.$$

(13)

with constants $\alpha ,\beta \in \mathbb{R}$.

The operational matrix (9) for this case is

$$\mathsf{D}={\mathsf{N}}^4-\alpha {\mathsf{M}}^2{\mathsf{N}}^2-2\alpha \mathsf{M}\mathsf{N}+\beta {\mathsf{M}}^4-(\alpha +\lambda )I$$

and ${\mathsf{B}}_n=[v_0;v_0{\mathsf{N}}^2;v_1;v_1{\mathsf{N}}^2]$, with the same $v_0$ and $v_1$ vectors of the previous example.

If ${\lambda }_{n,k}$ is the kth root of $det({\mathsf{T}}_n)$, and considering ${\tilde{\delta }}_{n,k}=\frac{{\lambda }_{n,k}-{\lambda }_{n-1,k}}{{\lambda }_{n,k}}$ as an estimative of the relative error in ${\lambda }_{n-1,k}$, then ${\delta }_n={max}_{k=1,\dots ,m}|{\tilde{\delta }}_{n,k}|$ is an estimative of the maximum relative error in the first m eigenvalues of the problem. In Figure 3 left, we present ${\delta }_n$, with $m=6$ and with $m=8$ for Example 3 with $\alpha =0.02$ and $\beta =0.0001$ for $n=16,\dots ,45$. In Figure 3 right, the absolute relative error $|{\tilde{\delta }}_{n,1}|$ in the lowest eigenvalue is presented for the same n values.

In

Table 2. Eigenvalues of Example 3 presented in [8] and [2] and ${\lambda }_{n,k}$ with $n=35$.

k	${\mathit{\lambda }}_{\mathit{k}}$ [8]	${\mathit{\lambda }}_{\mathit{k}}$ [2]	${\mathit{\lambda }}_{\mathit{n},\mathit{k}}$
1	0.21505086437	0.21505086436971492	0.2150508643697136
2	2.75480993468	2.754809934682884	2.754809934683077
3	13.2153515406	13.215351540558812	13.21535154055782
4	40.9508193487	40.95081975913755	40.95081975916199
5		99.05347813813880	99.05347806349896
6		204.35449348957833	204.3557322681771

, we compare our results with those of [2] for the first six eigenvalues, and of [8] for the first 4, obtained with values $\alpha =0.02$ and $\beta =0.0001$.

Example 4.

Now, we consider the Orr–Sommerfeld problem

$$\left\{\begin{array}{c}{y}^{(4)}(x)-2{\alpha }^2{y}^{″}+{\alpha }^{4}y=i\alpha R[(U-\lambda )(y^{″}-{\alpha }^{2}y)-U^{″}y],-1<x<1,\hfill \\ y(\pm 1)=y^{\prime }(\pm 1)=0,\hfill \end{array}\right.$$

(14)

with fixed constants $\alpha ,R$ and function U.

The particular case $U=1-x^2$ is the Poiseuille flow and, with $\alpha =1$ and $R=10000$ was treated in [3,4,5,12]. The operational matrix in that case is

$$\mathsf{D}={\mathsf{N}}^4-[(2{\alpha }^2+(1-\lambda )i\alpha R)I-i\alpha R{\mathsf{M}}^2]{\mathsf{N}}^2+({\alpha }^4+(1-\lambda )i{\alpha }^{3}R-2i\alpha R)I-i{\alpha }^{3}R{\mathsf{M}}^2.$$

Like in Example 2, this is an upper Hessenberg matrix with zeros intercalating its non-zero elements and we can reduce the problem dimension, splitting in two matrices the Schur complement ${\mathsf{C}}_n$ of the resulting Tau matrix ${\mathsf{T}}_n$. Choosing Chebyshev basis, this is the operational version of the Tau procedure of [5], where the author was confined to eigenvalues associated with symmetric eigenfunctions, which is equivalent to finding the eigenvalues of ${\mathsf{C}}_E$.

In [5], the author obtained for ${\lambda }_1=0.23752649+0.00373967i$ as an 8 decimal places exact value for the most unstable mode of this problem. Working with double-precision arithmetic, we obtain ${\lambda }_1=0.2375264889204038+0.003739670740170985i$. This value results with $n=58$ that is an $29\times 29$ matrix ${\mathsf{C}}_E$, the same dimension used in [5].

In addition, with $R_c=5772.22$, the smallest value of R for which an unstable eigenmode exists [5], and ${\alpha }_c=1.02056$, we get the results presented in

Table 3. Values of ${\lambda }_1$ of Example 4 with critical values $R_c=5772.22$ and $\alpha =1.02056$.

n	${\lambda }_1$
50	$0.2640017404987603-3.099641201518763i\times {10}^{-9}$
60	$0.2640017396216782-3.035288544655118i\times {10}^{-9}$
80	$0.2640017390806805-2.028898296212456i\times {10}^{-9}$
n	${\lambda }_1$ [5]
44	$0.26400174-1.7i\times {10}^{-9}$
50	$0.26400174+5.9i\times {10}^{-10}$

, together with those of [5].

4. Non-Polynomial Coefficients

In the previous section, we solved problems in the conditions of Corollary 1, i.e., with differential operators acting in polynomial spaces. In a more general situation, if some of the coefficients $p_i$ in (6) are non-polynomial functions, then the corresponding matrices $p_i(\mathsf{M})$ are functions of $\mathsf{M}$ instead of polynomial expressions.

If a non-polynomial function $p_i$ in (6) can be defined implicitly by a differential problem, with polynomial coefficients, then we can first of all use the Tau method to find a polynomial approximation ${\tilde{p}}_i$ to $p_i$ and use ${\tilde{p}}_i(\mathsf{M})$ to approximate the matrix $p_i(\mathsf{M})$.

Example 5.

Mathieu’s equation appears related to wave equations in elliptic cylinders [19]. For an arbitrary parameter q, the problem is to find the values of λ for which non-trivial solutions of

$$y^{″}(x)+(\lambda -2qcos(2x))y(x)=0$$

(15)

exist with prescribed boundary conditions.

It can be shown that there exists a countably infinite set of eigenvalues $a_r$ associated with even periodic eigenfunctions and a countably infinite set of eigenvalues $b_r$ associated with odd periodic eigenfunctions [19]. We are interested in reproducing some of those values given in there.

The operational matrix for problem (15) is

$$\mathsf{D}={\mathsf{N}}^2+\lambda I-2qcos(2\mathsf{M}).$$

Our first step to approximate Mathieu’s eigenvalues is to approximate matrix $cos(2\mathsf{M})$. This can be done by, firstly, considering the function $z(x)=cos(2x)$ as the solution of a differential problem, using Tau method to get a polynomial approximation $z_n(x)\approx z(x)$. In a second step, the operational matrix $\mathsf{D}$ is approximated by

$$\tilde{\mathsf{D}}={\mathsf{N}}^2+\lambda I-2q{z}_n(\mathsf{M})$$

and, finally, the last step consists in building the Tau matrix ${\mathsf{T}}_n$ and evaluating the zeros of its determinant.

We take integer values $q=0,1,\dots ,16$ and boundary conditions $y^{\prime }(-1)=y^{\prime }(\pi /2)=0$ to get $a_r(q)$ for even r, $y^{\prime }(-1)=y(\pi /2)=0$ for odd r, and $y(-1)=y(\pi /2)=0$ to get $b_r(q)$ for odd r and $y(-1)=y^{\prime }(\pi /2)=0$ for even r.

In Figure 4, we show Mathieu eigenvalues $a_r(q),r=0,\dots ,5$ and $b_r(q),r=1,\dots ,6$. Values were obtained with a 18th degree polynomial approximation $z_{18}\approx cos(2x)$ and a $36\times 36$ Tau matrix ${\mathsf{T}}_{36}$ in Chebyshev polynomials.

We can observe, as pointed out in [19], that, for a fixed $q>0$, we have $a_0<b_1<a_1<b2<\cdots$ and that $a_r(q),b_r(q)$ approach $r^2$ as q approaches zero.

Example 6.

Mathieu’s equation also appears coupled with a modified Mathieu’s equation in systems of differential equations as multi parameter eigenvalues problems. The particular case

$$\left\{\begin{array}{c}{y}_1^{″}(x_1)+(\lambda -2qcos(2x_1))y_1(x_1)=0,0<x_1<\frac{\pi }{2},\hfill \\ y_1^{\prime }(0)=y_1^{\prime }(\frac{\pi }{2})=0,\hfill \\ y_2^{″}(x_2)-(\lambda -2qcosh(2x_2))y_2(x_2)=0,0<x_2<2,\hfill \\ y_2^{\prime }(0)=y_2(2)=0\hfill \end{array}\right.$$

(16)

is studied in [6,7] associated with the eigenfrequencies of an elliptic membrane with semi axes $\alpha =cosh(2)$ and $\beta =sinh(2)$.

To approximate eigenvalues for this problem, we first have to approximate $cos(2x)$ and $cosh(2x)$ by polynomials. Considering, as in the previous example, $z_{16}\approx cos(2x)$ the 16th degree Tau solution of

$$\left\{\begin{array}{c}{z}^{″}(x)+4z(x)=0,0<x<\frac{\pi }{2},\hfill \\ z(0)=1,z(\frac{\pi }{2})=-1,\hfill \end{array}\right.$$

(17)

and $w_{16}\approx cosh(2x)$ as the same degree Tau solution of

$$\left\{\begin{array}{c}{w}^{″}(x)-4w(x)=0,0<x<2,\hfill \\ w(0)=1,w^{\prime }(0)=0,\hfill \end{array}\right.$$

(18)

then

$${\tilde{\mathsf{D}}}_1={\mathsf{N}}^2+\lambda I-2q{z}_{16}(\mathsf{M})$$

and

$${\tilde{\mathsf{D}}}_2={\mathsf{N}}^2-\lambda I+2q{w}_{16}(\mathsf{M})$$

are matrices approximating the operational matrices associated with differential equations (16).

For each fixed q, we define matrices Tau $T_1$ and $T_2$, representing Mathieu and modified Mathieu equations, respectively. Defining $a_n(q)$ the nth eigenvalue of $T_1$, in ascending order, and $A_m(q)$ the mth eigenvalue of $T_2$, in descending order, in [6], it was proved that $a_n(q)$ and $A_m(q)$ are analytical functions of q. Moreover, for each pair $(m,n)$, it was proved the existence and uniqueness of an intersection point of curves $a_n(q)$ and $A_m(q)$. Those intersections identify the eigenmodes of the elliptic membrane.

In Figure 5, we recover, and extend, figures presented in [6] and in [7]. Intersection points of $a_n(q)$, the almost vertical curves, and $A_m(q)$, the oblique curves, are the eigenpairs $(a,q)$ of (16).

5. Nonlinear Eigenvalues Problem

In some differential problems, the eigenvalues can arise in a nonlinear relation with eigenfunctions. Let us consider the following second order problem, related to Weber’s equation:

Example 7.

$$\left\{\begin{array}{c}{y}^{″}(x)+(\lambda +{\lambda }^2{x}^2)y(x)=0,-1<x<1,\hfill \\ y(-1)=y(1)=0.\hfill \end{array}\right.$$

(19)

The operational matrix corresponding to the differential equation is

$$\mathsf{D}={\mathsf{N}}^2+\lambda I+{\lambda }^2{\mathsf{M}}^2$$

and so $det({\mathsf{T}}_n)$ is a polynomial with degree $2(n-2)$ in $\lambda$.

In

Table 4. Eigenvalues of Example 7 with $n=30$ and $n=31$. Decimal places presented are those that coincide, until to the first distinct two.

k	${\mathit{\lambda }}_{\mathit{k}}(\mathit{n}=30)$	${\mathit{\lambda }}_{\mathit{k}}(\mathit{n}=31)$
$-5$	$-19.674904478$	$-19.674904482$
$-4$	$-13.62505355977$	$-13.62505355969$
$-3$	$-13.200062264051$	$-13.200062264066$
$-2$	$-7.0356879747644$	$-7.0356879747642$
$-1$	$-6.59716200235713$	$-6.59716200235723$
1	1.951702364990329	1.951702364990324
2	4.28611106118016	4.28611106118021
3	7.5459203349991	7.5459203349981
4	10.126005915959	10.126005915966
5	13.52870217426	13.52870217408

, we present the 10 eigenvalues closest to zero, obtained with $n=30$ and with $n=31$. We can verify that ${max}_{k=-5,\dots ,5}|{\lambda }_k(30)-{\lambda }_k(31)|<5\times {10}^{-9}$.

6. Conclusions

Since the pioneering works of Orzag [5] and Ortiz and Samara [9], the Tau method has been scarcely used to solve differential eigenvalues’ problems. With our work, we conclude that the Tau method is a competitive one if we want to evaluate with high accuracy the first eigenvalues, in a large kind of differential problem.