\(^\ast \)Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy
Str. Fântânele 57, 400320 Cluj-Napoca, Romania, e-mail: ghcalin@ictp.acad.ro.

\(^\bullet \)School of Mathematics, University of the Witwatersrand, Johannesburg, South Africa, e-mail: Bertin.Zinsou@wits.ac.za.

1 Introduction

The aim of this paper is twofold. We investigate analytically as well as numerically an elliptic one-dimensional second order eigenvalue problem with interior transmission conditions. The problem is a Dirichlet one, self-adjoint, involving a discontinuous coefficient. Thus we first show that eigenvalues are real and geometrically simple and the eigenvectors are orthogonal.

Then, using the FEM with linear test and trial bases (the so called linear hat functions), we find out the whole spectrum of the problem i.e. the set of all eigenpairs (eigenvalues and eigenvectors). The method is local and produces some generalised eigenvalue problems (GEP) containing pencils of sparse (tridiagonal) and possibly symmetric matrices. We also investigate the accuracy of numerical eigenpairs. The accuracy of the first (smallest) eigenvalues is estimated by computing their relative drift with respect to the order of approximation as well as with respect to the physical parameter involved in the definition of transmission conditions. As we could not find in literature numerical result for the problem at hand, we have tried to back them by some results obtained using a higher order method. Thus we use a collocation method based in Chebyshev polynomials. This is a very efficient global method which works equally well when the transmission conditions are not taken into account. In order to show robustness of spectral collocation, we solve an additional periodic SL problem with divergence coefficient exhibiting two discontinuities. This time periodicity requires a spectral collocation based on classical Fourier system (FsC).

Whenever the coefficient in the divergence form of the problem is continuous, the spectral method fairly works, i.e., the eigenvalues are obtained close to the machine precision (32 digits) and the convergence in computing eigenvectors is super-geometric. Due to the discontinuity of divergence coefficient, this convergence reduces to an algebraic one.

Actually we observe that our numerical results are in accordance with the theoretical (analytic) ones. Moreover, we try to make them more precise, for instance, to observe to what extent the discontinuity of coefficients influences the smoothness of eigenmodes.

Differential equations of the elliptic type involving discontinuous coefficients arise in many fields of science and engineering (see for instance [ 6 ] , [ 10 ] , [ 13 ] , [ 21 ] and [ 22 ] to quote but a few). For example, they model the dependence of the electrostatic potential on the electric charge in a conductor. They also model static deformation in elastic solids. The coefficients appearing in the equations are material parameters. Real materials are often composite and the material parameters then exhibit jumps across internal surfaces. These surfaces are often called interfaces. When interfaces are present, solutions to the PDEs must satisfy extra conditions, called transmission conditions. The presence of jumps at interior interfaces can cause scattering of waves, such as in a vibrating string made of two different types of materials. Transmission problems arise often in optics as well as various physics and electrical engineering problems.

Studying PDEs with jumps in the coefficients are challenging and the corresponding numerical methods used to solve them are more complex, than in the case of continuous coefficients (see for instance [ 2 ] ). We solve in this paper a model transmission problem in one dimension on a finite interval, imposing Dirichlet boundary conditions at the ends of the interval, as well as a periodic Sturm-Liouville problem involving discontinuous coefficients.

We can trace back the direct sum spaces method to the old paper of Zettl [ 24 ] as well as to a more recent paper [ 9 ] . We are perfectly conscious that in spite of the fact that the divergence coefficient \(p\left(x\right)\) is discontinuous its inverse \(1/p\) is locally integrable and thus the classical theory of existence for Sturm-Liouville eigenproblem is applicable (see for instance the Pryce’s monograph [ 20 ] and that of Zettl [ 25 ] ). Using the above mentioned method, we want to open another perspective for the class of problems at hand.

In [ 1 ] (see also [ 17 ] ) Aydemir and Mukhtarov investigate various qualitative properties of eigenvalues and corresponding eigenfunctions of a Sturm-Liouville problem with an interior singular point. They introduce a special Hilbert space along with a Rayleigh-Ritz formulation. We will introduce in Section  3 a slightly more general and useful weak formulation than this Rayleigh-Ritz formulation. However, Yserentant in [ 26 ] shows that for such formulation the computed eigenvalues as well as the corresponding eigenvectors converge to their continuous counterparts.

The organisation of the paper is as follows: In Section  2, we find out the most important analytic results concerning the eigenvalue problem at hand. In Section  3, we provide the weak (variational) formulation of the second order elliptic eigenproblem with transmission conditions. In Section  4, we summarise the simplex FEM as well as ChC for our problem. In Section  5, we provide the extensive numerical results obtained, i.e., the set of the first six eigenvalues and of the first six eigenvectors. We also advocate on the accuracy of our outcomes and comment on the dependence of the spectrum on the magnitude of jump. In addition, we solve a periodic eigenproblem with discontinuous coefficients in Subsection  5.1 using FsC. Consequently we observe that both spectral schemes, ChC as well as FsC, can be some serious competitors for the well established codes devoted to the eigenvalue problems. Section  6, ends up the paper with some conclusions and open problems.

2 Eigenvalues and eigenfunctions

We consider the Sturm-Liouville (SL) equation,

where

With the differential equation 1, we consider the boundary conditions

together with the transmission conditions

The Hilbert spaces \(L_2([0, 1/2))\) and \(L_2((1/2,1])\) are the spaces of all classes of complex-valued measurable functions \(f\) and \(g\) such that

Define

and the inner product

for all \(f,g\in H\). Then \(H\) is a Hilbert space with respect to the inner product 9.

Since \(\tau (u)=\lambda u\) and \(\tau (v)=\mu v\), then it follows from the Lagrange identity and Green’s formula that

The boundary conditions 4 and 5 give \([u,v](0)=[u,v](1)=0\), while the transmission conditions yield \([u,v](1/2^-)=c^2[u,v](1/2^+)\). Hence \((\lambda -\overline\mu )(u,v)=0\); so if \(\lambda \ne \overline\mu \), then \((u,v)=0\).

Let \((\lambda _0,u_0)\) be any eigenpair of the problem 1–7. The pair \((\overline\lambda _0,\overline u_0)\) is an eigenpair of the complex-conjugate of the problem 1–7. Thus

and

Putting 11 and 12 into 10, we get \((\lambda _0-\overline\lambda _0)||u_0||^2=0\). Hence \(\lambda _0=\overline\lambda _0\) and the proof is complete.

Assume that there exist two linearly independent eigenfunctions \(u_0\) and \(v_0\) for the same eigenvalue \(\lambda _0\). The boundary condition 4 implies that \([u_0,v_0](0)=0\). Thus \([u_0,v_0](x)=0\) for all \(x\in [0,1/2)\). Since \(u_0\) and \(v_0\) are solutions of 1, then there exists \(\alpha _1\) such that \(u_0(x)=\alpha _1v_0(x)\) for all \(x\in [0,1/2)\). Similarly, it follows from the boundary condition 4, that there exists \(\alpha _2\ne 0\) such that \(u_0(x)=\alpha _2v_0(x)\) for all \(x\in (1/2, 1]\). Hence

Substituting 13 in the transmission conditions 6–7, we get

It follows that \(\alpha _1-\alpha _2=0\). Therefore \(u_0\) and \(v_0\) are linearly dependent on \([0,1/2)\cup (1/2,1]\). This completes the proof.

Lemma 1, Theorem 2 and Remark 5 lead to

3 Variational (weak) formulation

From the definition of coefficient \(\alpha \) (see Section 2) there exists a positive constant \(\gamma \) such that \(\alpha \left( x\right) {\gt}\gamma {\gt}0,\) i.e., the problem is elliptic. To justify the transmission conditions, we begin with the so-called weak formulation of the problem (1)–(5), namely:

Such formulations for elliptic problems are the starting point for the most modern numerical methods (see the classic texts of Brenner and Scott [ 7 ] or Ciarlet [ 8 ] ).

In order to obtain the equation (18), in the unknown \(u\), we multiply the equation (1) with another suitable function \(v\), called a test function, and then integrate over \(\left[0,1\right]\). Thus we simply have

The function \(v\) is chosen to satisfy the same Dirichlet boundary conditions as \(u\), namely:

so that there are no extra terms coming from the integration by parts.

Then, we split up the integral on the left hand side of (19), given that \(\left( \alpha \left( x\right) u^{\prime }\left( x\right) \right) ^{\prime }\) has a jump at \(x = 1/2\) and get

This allows us to have continuous functions in each integrand. We can then integrate by parts each integral, obtaining the following:

We can see given the boundary conditions that \(\left[ u^{\prime }v\right] _{x=0}\) and similarly \(\left[ u^{\prime }v\right] ^{x=1}\) will vanish. In order to obtain the weak formulation (18), we need to impose some conditions at the jump site \(x=1/2\). More exactly we must apply first the continuity of the solution \(u\) at this site. Thus, we need to enforce the transmission condition (6). Hence we have to impose just the transmission condition for the first derivative (7). Thus, we obtain (18) which is the weak formulation of the transmission eigenproblem (1)–(7). However, in this formulation the coefficient \(\alpha \) is not differentiated and the integral is well defined as long as the test (shape) function \(u\) and the trial function \(v\) are differentiable in some sense.

The formal variational (weak) formulation (18) now reads:

Find \(u\in H_{0}^{1}\left(0,1\right) \) and \(\lambda \in \mathbb {R}\) such that

This weak solution \(u\) with some supplementary assumptions is also the strong (classical) solution to (1)–(7) (see for instance the classical monographs of Brenner and Scott [ 7 ] , Ciarlet [ 8 ] or Nec̆as [ 18 ] for rigorous aspects of variational calculus in Sobolev spaces).

4 The Simplex FEM and ChC

Let now \(\mathcal{P}_{1}\) be the space of polynomials of degree at most \(1\). We define the space of piecewise linear polynomial functions

We will use the FEM in its simplest form, i.e., simplex where the basis for test as well as trial functions, contains only piecewise linear functions, the so called linear hat functions. All these functions belong to the above defined space \(X_{1}\).

In order to define these functions we start by partitioning the interval \([0, 1]\) into \(N\) equal subintervals, where \(N\) is a positive integer called the approximation order. This gives rise to \(N+1\) nodes of the form \(x_{j}:=j/N, j=0,1,\ldots ,N\). Later on, we will choose \(N := 2k\) for some integer \(k\), so that \(1/2\) is a node and the jump in the coefficient \(\alpha \) occurs at that node. This simplifies our analysis as we discuss more later on.

For each \(1\leq j \leq N-1\) we construct a hat function \(\phi _{j}\) as follows. Each \(\phi _{j}\) is a hat function of height \(1\), which is non zero and piecewise linear on the interval \(\frac{j-1}{N} {\lt} x{\lt} \frac{j+1}{N}\). Thus,

and its derivative reads

We will seek a solution to our problem (20) in the form:

The stiffness matrix \(\mathbf{A}\) has the entries

and thus it is a tridiagonal one.

The mass matrix \(\mathbf{M}\) is defined by the entries

This matrix is also tridiagonal and additionally symmetric. Thus, the non zero entries in the stiffness matrix \(\mathbf{A}\) and in the mass matrix \(\mathbf{M}\) have the following numerical values:

and

We plug in the computed entries \(A_{i,j}\) and \(M_{i,j}\) to get the following sparse GEP

Being sparse and tridiagonal, both matrices in GEP (26) can be efficiently implemented in MATLAB using the routine diag.

The analytic results from Sect.  2 and Sect.  3 concerning the continuity of \(u\), i.e., \(u\in H_{0}^{1}\left(0,1\right)\) enable us to use alternatively the strong ChC method in order to solve the SL transmission problem. Thus in the finite dimensional representation of solution \(u\) in (??) instead of linear hat functions we use the Chebyshev polynomials (see for instance our text [ 12 ] and the seminal paper [ 23 ] ).

Thus, we have to find out the eigenpairs of the matrix

the vector \( \mathbf{X} \) contains Chebyshev nodes of the second kind \(x_{k}, k=1,\ldots ,N\) and the \(N\) dimensional vector \(\alpha \left( \mathbf{X}\right)\) is defined as

The matrix \(\mathbf{D}_{ChC}\) is the Chebyshev collocation differentiation matrix on the above Chebyshev nodes (see [ 23 ] for its implementation). Actually we search the eigenvalues of

as we have enforced the homogeneous Dirichlet boundary conditions. The MATLAB code eig and another real variant of the Jacobi-Davidson method (see our contribution [ 11 ] ) are used consecutively in order to mutually confirm the numerical values obtained. We have to mention that the matrix

is fully populated and rather non normal. Its Henrici’s number equals

For the importance of non normality and its numerical measure provided by Henrici’s number in context of eigenvalue problems we refer to our
contribution  [ 12 ] .

5 Numerical results

A banded FEM (Galerkin) discretization matrix saves us nothing over a dense matrix in pseudospectral method if the linear algebra is to be handed off by QZ (QR) algorithm. Its main drawback is the cost. Because it uses iteration, rather than a finite set of steps, a precise estimate of cost is impossible but experience has shown that the QZ cost is \(O(10 {N^3})\) floating point operations. Alternatively, we use the MATLAB built in eigs which takes into account the sparsity of matrices in (26) and finds out the left column in Table 1.

For a specified eigenvalue \(j\) the relative drift is defined by J. P. Boyd in [ 5 ] with the quotient

where \(N_{1}\) and \(N_{2}\) are two distinct orders of approximation. The relative drift for the first \(25\) eigenvalues computed by FEM and by ChC is displayed in Fig. 1 panels A) and respectively B).

It signifies a reasonable numerical stability of computing eigenvalue process in FEM case and an excellent stability in case of ChC.

Fairly interesting is also the relative drift of eigenvalues with respect to the exact ones, i.e., when \(c:=1\). In this case, we know that the exact eigenvalues of (1)–(5) are \((\pi k)^2\), \(k=1,2,...\). The drift in this case is depicted in Fig. 2. It is clear from panel B) of Fig. 2 that ChC computes the first few eigenvalues close to the machine precision.

The accuracy of FEM is clearly much worse, i.e., of \(O(h)\) (see panels A) of Fig. 1 and 2).

Moreover, using a multi-precision computing toolbox [ 15 ] , our numerical experiments show that the accuracy of ChC increases with the number of digits in multi-precision approximation. As it is apparent from panel A) of Fig. 2 for FEM, this does not happen.

Actually, working with \(100\) digits precision the first eigenvalue is computed by FEM with an accuracy better that \(O(10^{-80})\).

It is also interesting to observe the asymptotic behaviour of the computed eigenvalues. In the Zettl’s monograph [ 25 , p. 73 ] , it is known that

When we try to verify this statement we get Fig. 3.

This picture confirms the Boyd’s RULE-of-THUMB for eigenvalues, i.e., for a spectral method using \(N + 1\) terms the lowest \(N/2\) eigenvalues are usually accurate within a few percent while the larger \(N/2\) are useless (see [ 5 ] , p.132).

Fortunately, but less expectantly, a fairly similar conclusion holds for the set of eigenvalues computed by FEM.

The first six eigenmodes computed by ChC are displayed in Fig. 4. We observe that they are continuous at the transmission point but their derivatives are not.

In order to evaluate the asymptotic order of convergence of ChC method in computing eigenmodes, we use the fast Chebyshev transform. Thus, we get the coefficients of Chebyshev expansion and then plot their absolute values in a log-linear plot. The behaviour of these coefficients for the first four vectors is depicted in Fig. 5. Using the strategy from Boyd [ 5 ] we can imagine the so called envelope, i.e., a curve which bounds these coefficient from above. This indicates an algebraic convergence, i.e., the coefficients \(\widehat{u}_{k}\) satisfy \(\widehat{u}_{k}\sim k^{-2}\) for large \(k\).

As we compute a large set of eigenfunctions, Corollary 6 from Sect. 2 provides us an excellent tool to assess the accuracy in their computing. Thus, in order to measure the departure from orthogonality of the computed eigenmodes in a log-linear plot, we depict the numerical values of discrete scalar product of eigenvectors, namely

5.1 SL problem with discontinuous coefficients

In order to underline the versatility of spectral collocation methods we briefly compute the eigenspectrum of a SL problem with discontinuous coefficients from the paper of Babus̆ka and Osborn [ 2 ] . They consider a \(2\pi \) periodic boundary value problem attached to equation (1) when the divergence coefficient \(\alpha \) is measurable. The problem is factorised and a variational formulation is used. Some convergence results and error estimates for a mixed FEM are derived. These error estimates are based on the application of Sobolev spaces with variable constant order. Unfortunately no numerical results are provided.

However, the coefficient \(\alpha \) is defined by

We solve this periodic problem by Fourier collocation. It means that we have to use for the discretization the first order Fourier collocation differentiation matrix on the equispaced nodes

(see again [ 23 ] Sect. 3.5). This differentiation matrix is now a circulant one and thus normal. The first six eigenvalues are displayed in Table 2.

We did not find similar results in the literature in order to validate the ones in the Table 2.

In order to evaluate the asymptotic order of convergence of FC method in computing eigenmodes, we use the discrete FFT (fft from MATLAB). Thus, we find out the coefficients of Fourier expansions. A fairly similar analysis with that for ChC method shows the same order of convergence, i.e. algebraic (see lower panels in Fig. 7).

It is worth noting at this moment that in [ 2 ] the authors solve only theoretically this problem and prove that the rate of convergence in approximation by FEM based on trigonometric polynomials is of order \(N^{-2}\) and this estimate cannot, in general, be improved. Our numerical outcomes reported in the lower panels of Fig. 7 confirm this statement.

Eventually, we have to observe that the sets of numerical eigenvalues computed by FEM, as well as spectral collocation methods, satisfy by far the estimations for lower bounds provided in [ 14 ] .

6 Brief conclusions and some open problems

On completely different analytical considerations than the classical ones, we have proved three essential properties for the SL transmission eigenproblem, namely: all the eigenvalues of the problem are real, each eigenvalue is geometrically simple and eigenfunctions corresponding to distinct eigenvalues are orthogonal. To solve the problem numerically we have resorted to two distinct methods. The results obtained with their help are reasonably mutually closed. Additionally we remark that the global ChC method is by far much more precise than FEM despite the fact that it is based on full populated and non normal matrices compared to the sparse (even symmetric) matrices produced by the FEM.

It is also important to underline that the global collocation method operates without considering the transmission conditions. It is robust, efficient and easy to extend to higher-order problems.

The Fourier collocation worked with the same convergence order for a periodic problem. Last but not least important, both spectral collocation methods offer challenging alternatives for older codes devoted to solve various SL problems (see for instance [ 3 ] , [ 4 ] and [ 19 ] ). Finally, we must also note the reciprocal confirmation of the analytical and numerical outcomes.