Numerical Approaches to Compute Spectra of Non-Self Adjoint Operators and Quadratic Pencils

Fatima Aboud,Fran?ois Jauberteau,Guy Moebs and Didier Robert

1 Mathematics Department,College of Science,University of Diyala,Baquba,Iraq

2 Laboratoire de Mathématiques Jean Leray,CNRS-UMR 6629,Université de Nantes,France.

Abstract.In this article we are interested in the numerical computation of spectra of non-self adjoint quadratic operators.This leads to solve nonlinear eigenvalue problems.We begin with a review of theoretical results for the spectra of quadratic operators,especially for the Schr?dinger pencils.Then we present the numerical methods developed to compute the spectra:spectral methods and finite difference discretization,in infinite or in bounded domains.The numerical results obtained are analyzed and compared with the theoretical results.The main difficulty here is that we have to compute eigenvalues of strongly non-self-adjoint operators which are very unstable.

Key words:Nonlinear eigenvalue problems,spectra,pseudospectra,finite difference methods,Galerkin spectral method,Hermite functions.

1 Introduction

We are interested here in equations likeL(λ)u=0 whereL(λ)is a linear operator on some linear spaceE,depending on a complex parameterλ.WhenL(λ)=L0?λI,this is the usual eigenvalue problem:findλ ∈C andu ∈E,0 such thatL(λ)u=0. In many applications,in particular for dissipative problems in mechanics,it is necessary to consider more general dependence in the complex parameterλ.A particular interesting case is a quadratic family of operators:L(λ)=λ2L2+λL1+L0,whereL2,L1andL0are linear operators inE.We shall say thatL(λ)is a quadratic pencil.

Let us consider the second order differential equation:

Eq. (1.1)is a model in mechanics for small oscillations of a continuum system in the presence of an impedance force[26].Now by looking for stationary solutions of(1.1),u(t)=u0eλt,we have the following equation:

So Eq.(1.2)is a non-linear eigenvalue problem in the spectral parameterλ ∈C.We say thatλis a non-linear eigenvalue if there existsu00 satisfying(1.2).

The operatorL1represents a damping term as we see in the following simple example.Let us consider the perturbed wave equation(see[5]):

wheret∈R andx∈T:=R/2πZ.Here we haveL2=IE(identity operator),L1=?2aandL0=??2/?x2.The damping terma ＜0 is here constant.So we have to solve(1.3)with periodical boundary conditions.The stationary problem is reduced to the equation:

Then we have fork2≥a2the damped solutions of(1.3):

Whenais a function ofxwe have no explicit formula so we need numerical approximations to compute the damping modes.It is the main goal of this work.

Such generalized eigenvalue problems have appeared in a completely different way.The question was to decide if a class of P.D.E with analytic coefficients preserves or not the analyticity property.To be more explicit,let us consider a P.D.E:Pu=f.Assume thatfis analytic in some open set Ω,is-it true thatuis analytic in Ω?This is true for elliptic operators.For some example,this question can be reduced to the following(see[23]for more details):does there exist(R)such that

Existence of non null solutions for(1.2)and(1.4)is a non trivial problem.For(1.4)it was solved in[33]where it is proved that the generalized eigenfunctions span the Hilbert spaceL2(R).

On the other side we can prove that the equation:

has only the trivial solution0 inL2(R),?λ∈C.

Our aim in this work is to present several numerical approaches concerning this kind of non-linear eigenvalue problems.For simplicity we only consider quadratic pencil such thatL0=I.We can reduce to this case ifL0orL2are invertible in the linear spaceE.To every quadratic pencilL(λ)we can associate a linear operatorALinE×Esuch thatλis a non-linear eigenvalue forLif and only ifλis a usual eigenvalue forAL.The operatorALis called a linearization ofL(λ).It is easy to see that we can choose:

So non-linear eigenvalue problems(for polynomial operator pencils)can be reduced to usual eigenvalue problems but it is useful to take care of their particular structure.There exists infinitely many linearizations.

We are mainly interested here in the case called Schr?dinger pencils:

whereVandaare real smooth functions from IRminto IR(see Section 2.2 for more details).The main questions we want to discuss is the localization of the eigenvalues ofLV,ain the complex plane C.

In Section 3 we discuss several numerical approaches for the computation of the spectra of quadratic pencil,in infinite and in bounded domain.Using spectral methods or finite difference discretization to approximate the operatorsL0,L1andL2,we obtain a nonlinear eigenvalue problem.Then,after linearization we obtain an eigenvalue problem to solve.The numerical results obtained are analyzed and comparisons with theoretical results are done.The main difficulty is that we have to compute eigenvalues of strongly non self-adjoint operators,which leads after discretization to strongly unstable nonlinear eigenvalue problems.Finally,in Section 4 we give conclusions and future works.

2 A review of results obtained by functional analysis methods

In this paper our main goal is to locate,in the complex plane,the generalized spectrum of pencils of differential operators.These generalized eigenvalues are very unstable and the theorem of existence of eigenvalues may be difficult to prove.Moreover for numerical results it is difficult to make the difference between the eigenvalues and pseudoeigenvalues.

Of course there exists a huge number of papers concerning computations of eigenvalues.Here we only notice some references directly connected with our paper.

For estimation of the spectra of linear differential operators we can refer to[3,7,8,11,12,22,28,32].For numerical methods to solve nonlinear eigenvalue problems(finite dimensional case)and for linearization of such nonlinear eigenvalue problems we can refer to[5,6,9,19,21,31,36,39].

2.1 Quadratic family operators

For more details on the results presented in this section we refer to the book[30].

Let us consider the quadratic family of operatorsL(λ)=L0+λL1+λ2whereL0,L1are operators in an Hilbert spaceH.IfHif of dimensionN ＜+∞the eigenvalues are the solutions of the polynomial equation det(L(λ))=0.WhenNis large this could be a difficult problem at least for numerical computations.In applications involving PDE,His aL2space or a Sobolev space,which is infinite dimensional and there is no explicit equation for the generalized eigenvalues.Moreover,as we shall see later,the non linear eigenvalue problem is equivalent to a linear eigenvalue problem which,in general,is non self-adjoint hence unstable.

Theorem 2.1.L(λ)is a family of closed operators in is meromorphic in the complex plane.The poles λj of L?1(λ),with multiplicity mj,coincide with the eigenvalues withthe same multiplicities,of the matrix operator AL in the Hilbert space ,with domainwhereLV,a(λ)satisfies the above theorem for

IfL0is positive and non degenerate we have the symmetric linearization:

Let us denote Sp[L]the eigenvalues ofAL(which coincide with the poles ofL?1(z)).

Ifλ0∈Sp[L]we denote byEL(λ0)the linear space of the solutions{u0,u1,···,uk,···}of the equations:

The dimension ofEL(λ0)is the multiplicity ofλ0(for details see[33]).

Assume thatL0,L1are self-adjoint,L0is positive non degenerate and that there existκ≥0 andδ≥0 such thatis a bounded operator onHandAssume thatis in the Schatten classCp,p≥1.

Theorem 2.2.If0＜δ≤1/2then the spectra of L is the domain:

=0dense in H.

Forδ＞0 we get that the eigenvalues are localized in a vertical parabolic domain in the imaginary direction.Forδ=0 endκsmall the eigenvalues are localized in a small sector around the imaginary axis.Notice that forκof order 1 the above theorem does not give any information on the location of Sp[L];we only know that it is a discrete and infinite subset of C.

Moreover ifL1has a sign we have:

Proposition 2.1.If L1≥0thenSp[L]?{λ∈C,?λ≤0}.If L1≤0thenSp[L]?{λ∈C,?λ≥0}.

The above result applies for example to:

For this example we havehence the spectra is localized inside the parabolic region

For Schr?dinger pencilsLV,awe can say more.

2.2 More results for Schr?dinger pencils

Let us recall our definition of Schr?dinger pencilsLV,a(λ)=?△+V?2aλ+λ2.

We assume that the pair of functions(V,a)satisfies the following technical conditions that we refer as[cond(V,a)](we do not try here to discuss the optimality of these conditions):

V,aare smoothC∞functions on Rdand there existsk＞0 such that

Under these conditions we know thatL0=?△+Vis an unbounded self-adjoint operator inL2(Rd)and for everyλ∈CLV,a(λ)is a closed and Fredholm operator with the following weighted Sobolev space domain:

Moreover the set Sp[LV,a]of eigenvalues ofLV,ais a discrete set(empty or not),each eigenvalue having a finite multiplicity and the only possible accumulation point in the complex plane is ∞.

Notice that ifλis an eigenvalue then its complex conjugateis also an eigenvalue.We have the following result:

Proposition 2.2.Assume that(V,a)satisfies[cond(V,a)]and that a ≤0,a(x0)＜0for some x0∈Rd.ThenSp[LV,a]is in the open sector

Proof.Letu ∈L2(Rd),such thatLV,a(λ)u=0.Setλ=r+is.We know thatr ≥0.Assume thatr=0.Reasoning by contradiction we first prove thats=0.Ifwe get thatvanishes in an non empty open set of Rdand applying the uniqueness Calderon theorem for second order elliptic equation we getu=0 on Rdand a contradiction.Ifs=0 we get:

So again we get thatuvanishes on a non empty open set and a contradiction like above.

Let us remark that the general results given in Theorem 2.2 apply if there existsδ≥0 such that|a|(x)≤CV(x)1/2?δor|a(x)|≤κV(x)1/2withκsmall enough.

For 1D Sch?dinger pencils accurate results were obtained by M.Christ[14,15]and by[13].Let us recall here some of their results.They consider the pencils

withk∈N.Here we shall only considerkeven.The above assumptions are satisfied.We have:

Theorem 2.3([13],Theorem 1).Sp[Lk]such that|λ1|＜···＜|λn|＜|λn+1|＜···.Then we have for n→+∞:

Theorem 2.4([33],[1]).The linear space spanned by the generalized eigenfunctions associated with the eigenvalues{λn}is dense in L2(R).

In[33]the proof was given forL2(λ)and forLk(λ),k＞2 andkeven in[1],[2].

In the following result we shall see that the spectral set Sp(Lk)is very unstable under perturbations.M.Christ[15]has considered the following model:

We also have

wherePis a polynomial.Assume that the degreekofPis even,+a1x+a0.We have:

Proposition 2.4.SpIn other words for every λ∈C??it is known that every solution in is in the Schwartz space S(R)(see[33])in the Schwartz space S(R).

3 Spectra approximation with spectral methods and FDM

The aim of this section is to present different numerical methods to analyze the localization in the complex plane of generalized spectrum of operator pencils like(2.1).More precisely we compare our numerical results with the theoretical results recalled in Section 2.Proposition 2.6 gives two sectors where one can find eigenvalues;Theorem 2.7 gives more accurate results on the localization for the eigenvalues with large modulus(when the modulus increases the eigenvalues are near the imaginary axis).The theoretical results of Section 2 are given in an infinite domain.But for numerical computations it is often needed to truncate the domain(for example for the finite difference method)for working in a finite domain.Then the question is to optimize the size of the domain to reduce the error term.This may be a delicate problem.Indeed the size of the domain can change dramatically the spectrum of the considered operator.For example takingk=1 in(2.1)we can see that the spectrum is empty in the infinite domain R.But considering any finite interval and homogeneous Dirichlet conditions one can prove that the spectrum of(2.1)is discrete and infinite.

For the kind of pencils of differential operators considered in this paper the eigenfunctions are fast decreasing at infinity. But an eigenfunctionφλassociated with the complex eigenvalueλis small only if|x|≥f(?λ)for some increasing functionfsuch that limr→∞f(r)=∞.So we have to consider larger and larger domain to capture eigenvalues with large imaginary part.

In a first step we are staying with an infinite domain and we develop a spectral method using Hermite functions.This method is very accurate for our problem because the eigenfunctions are in the Schwartz spaceS(R).

In a second step we consider the finite domain approximation(with homogeneous Dirichlet conditions)and apply a finite difference method or a spectral Legendre-Galerkin method.As already said we have to optimize the size of the domain.This is solved by using properties of the zeros of Hermite functions.

Notice that the spectrum of the considered operator pencils is the spectrum of nonself-adjoint operators which are far from self-adjoint or normal operators.So the spectrum is very instable and it seems useful and important to compare different numerical methods.We also compare for our models the numerical computations of the spectrum with computations of the pseudo-spectrum.

Hereafter we present numerical methods for the computation of the spectra of linear operators with quadratic dependence of the spectral parameterλ(quadratic pencil),see(1.2):

whereL0andL1are operators on some Hilbert spaceHand I is the identity operator.

3.1 A first quadratic pencil

In this section we consider the following operator:

whereais a real parameter.We want to solve the following nonlinear eigenvalue problem:

Fora=1 we recover the problem(1.4).Notice that ifusatisfies(3.2)then we know thatu∈S(R)(see[33]).

The problem(3.2)can be reformulated as an eigenvalue problem after linearization ofLa(λ)(see Section 1).Indeed,if we setv=λuwe can rewrite(3.2)as:

where

3.1.1 Eigenvalue computations with Hermite spectral method(unbounded domain)

We look for an approximationuNofusuch thatLa(λ)uN=0,with:

In order to computeuN,we use a method of weighted residuals(MWR,see,e.g.,([10,20]):

where(.,.)is the scalar product inL2(R). SettingvN=λuN,using the orthogonality properties of the Hermite function inL2(R)and the relations(A.2),we obtain after linearization ofL(λ)the following eigenvalue problem:

which is an approximation of the eigenvalue problem(3.3).UN(resp.VN)is the vector containing the coefficients(resp.ofu(resp.v),k=0,...,N.The matrixAa,Nis the square matrix of order(2N+2):

whereand=0,...,N,with

Note thatL0Nis a pentadiagonal symmetric matrix such that

forj=0,...,N,where

Moreover,is a tridiagonal symmetric matrix such that

For the numerical computation of the spectrum ofAa,Nwe use the function DGEEV of the LAPack library which is based on the Schur factorization.

Fora=1,in order to analyze the spectrum of the continuous operator(1.4),we consider a simplified operator,deduced from the operator(3.1)fora=1,wherexis replaced with a real constantb.We obtain the following problem:

We look for a solutionu(x)of the problem(3.6)of the form=0,...,N.Substituting in(3.6)and using the relations(A.2),we obtain:

Using the scalar product inL2(R)of(3.7)withφkgives

Figure 1:Spectrum of the matrix Aa,N(3.5)(Hermite spectral method)for N=50 and a=1.

So the size of the containment area is

and we can retainb=Las value forbin the estimation(3.9).Here,forN=50 following(3.10)we deduce that the size of the containment domain is 2Lwith.In Figure 1,we can see that effectively we have 0≤λN,r ≤L2.

3.1.2 Eigenvalue computations with finite difference method(bounded domain)

The operatorLa(λ)(see(3.1))is defined on the domainD(A)={u∈H2(IR),x4u∈L2(IR)}.Souis decreasing asO(1/x4)when|x|is increasing(see(A.6))and we want to consider the nonlinear eigenvalue problem(3.2)in a bounded domain with homogeneous Dirichlet boundary conditions:fnidsuch that

As before,the problem(3.11)can be reformulated as an eigenvalue problem

where

To obtain an approximation of the continuous problem(3.12)in a finite dimensional space,we consider on the domain Ω a mesh-grid with a mesh Δx=2L/Non Ω and we notexj=?L+jΔx,j=0,...Nthe points of the grid.We have retained homogeneous Dirichlet boundary conditions forx=±L,sou(x0)=u(xN)=0.We look for an approximationuN,vNofuandv=λusuch that:

withUNandVNtwo vectors containing respectively the approximationsuN(xj),vN(xj)ofu(xj),v(xj)andAa,Nis the square matrix of order 2N?2:

where

is the discretization of the operatorL0with a second order centered finite difference scheme and

L0Nis a tridiagonal symmetric matrix such that

andis a diagonal matrix such that

For the numerical computation of the spectrum of the matrixAa,Nwe use the function DGEEV of the LAPack library.

Now,we are interested to analyze the dependence of the spectrum of the operator(3.1)in function of the real parametera.For this,we consider an approximation of the infinite dimensional domain as a bounded domain with periodic boundary conditions.We look for eigenfunctionuk(x)=?kexp(ik′x),withk′=kπ/L,of the continuous operator(3.1).ComputingLa(λ)uk(x)we obtain the following equation:

Whenais increased from 0 to 1 the ratio of the imaginary part over the real part ofλ,

is decreased and it is infinite fora=0.We can observed this on the numerical simulations corresponding toa=0,a=0.5,a=0.9 anda=1.0 obtained with the finite difference scheme forN=50 andL=10(see Figure 2).

We are now interested with the operatorLa(λ)fora=1.We look for the spectrum of the discretized operator,using finite difference method,wherexis replaced with a real constantbas in(3.6). Indeed,when we discretize with a finite difference scheme,we consider thatxis constant over one spatial step Δx=2L/N.So,in a first approximation,we consider an operator deduced fromLa(λ)in whichxis chosen constant equal tobover all the domain Ω.We have:

If we consider periodic boundary conditions,we look for a solution of(3.14)of the formwithk′=kπ/L.Substituting in(3.14)and supposing thatwe obtain

Figure 2:Spectrum obtained with the finite difference scheme for a=0,a=0.5,a=0.9,a=1.0,N=50 and L=10.

Finally we haveλN=λN,r+iλN,iwithλN,r=b2is wavenumber independent and

is wavenumber dependent.So the spectrum of the discretized operator is located in the part of the complex plane such that

Figure 3:Spectrum of the matrix Aa,N(3.13)(finite difference scheme)for N=50,L=10 and a=1.

in agreement with|λ|≤L2.

Numerical instability

In order to study the numerical instability of the finite difference scheme in function of the mesh-gridxj,j=0,...N,we consider a small perturbation on each point of the grid,xj+ε,j=0,...N,whereεis a small parameter.The matrixAa,N(see(3.13))is replaced with the matrix:

whereENis a perturbation matrix.If we compare the eigenvaluesλNof the matrixAa,Nwith the eigenvaluesλN,εofAa,N,εwe have:

Figure 4:Condition number of the eigenvalues λN in function of the modulus|λN|for N=50,L=10 and a=1.

whereUN,εis a right eigenvector ofAa,N,ε.So we deduce(see[35]):

Forε=0 we obtain:

If we multiply on the left the previous equality withVNa left eigenvector ofAa,Nwe obtain:

whereandUNis a right eigenvector ofAa,N.The equality(3.15)measures the sensitivity of the eigenvalueλNof the matrixAa,Nin function of a perturbationεon the mesh-grid(condition number of the eigenvalueλN).

Here the matrixENis the matrix of order 2N?2:

whereE0,N(resp.E1,N)is the diagonal matrix with the elements(resp.4axj)on the diagonal,j=1,...,N?1(we have neglected inENthe terms inεn,withn＞1).

Figure 5:Average on the eigenvalues computed with the finite difference scheme using 11 staggered grids for a=1,N=1000 and L=10.

Figure 4 presents the condition number of the eigenvaluesλNin function of the modulus|λN|forN=50,L=10 anda=1.We can see that eigenvalues are ill conditioned,excepted for the eigenvalues with small modulus.This can explain the convergence problem whenNis increased.If we compare with the rotated harmonic oscillator the condition numbers of the eigenvalues are much greater for the operator(3.1)(around 108)than for the rotated harmonic oscillator(around 10).A small perturbation on the grid points induces large perturbations on the eigenvalue computations.

For the use of staggered meshes to avoid spectral pollution,we may mentioned the following reference[29].

Pseudospectra

Now we consider the pseudospectra since it is known that the numerical computation of the pseudospectra is more stable than for the spectra(see Appendix B).To obtain the pseudospectra,following Defniition B.2 we look forsuch thatis large,i.e.:

Figure 6:Computation of the pseudospectra(3.16)of the matrix Aa,N(3.13)(finite difference scheme),for N=1000,L=10 and a=1. The different colors correspond to different size of the small parameter ε. The smallest parameter ε corresponds to orange snowflakes.

When the mesh-grid retained in the complex plane is fine,the pseudospectra computation is expensive since it requires to compute the minimal singular valueat each pointzof the mesh.So we use parallel computation in order to accelerate the computation.The numerical solution is done thanks to the linear algebra library LAPack which contains specialized algorithms for singular values problems,especially the one called ZGESVD for complex matrices in double precision(ZGESVD is based on bidiagonal QR iteration).As the matrix(3.13)is quite huge,and computing time a bit long,a parallelization by MPI(Message Passing Interface)is implemented with the client/server model.One process(the server)distributes values of the complex parameterz(see(3.16))to the other processes(the clients)which sample the domain.The server renews their data as the work progresses. Each client builds the matrix to be study and sends to the server,at the end of the computation,the smallest value.This system has the advantage of being dynamically balanced.As there is no communication(in MPI sense)between the clients,the efficiency of the parallelization is complete. As an example,the simulation corresponding to the parametersN=5000,L=1000,a=1 and to an area of the complex plane[0,150]×[0,150]with a mesh stepdx=1 anddy=1 in the real and imaginary directions has needed 40 cores(Intel Xeon E5-2670 at 2.5GHz)during quite 40 days.

3.1.3 Eigenvalue computations with Legendre spectral Galerkin method

In order to obtain a higher accurate numerical scheme in bounded domain,we propose a spectral numerical scheme using Legendre Galerkin basis.

We consider the problem(3.11).This problem is reformulated as an eigenvalue problem(3.12).But instead of using a finite difference scheme to obtain an approximationuN,vNofuandv=λu,we use a spectral method with Legendre Galerkin basis Φl.Such basis is obtained as a linear combination of Legendre polynomials:

and

Moreover,we need the expressions ofx2Φlandx4Φlas linear combination of the Legendre polynomials.We have:

with

In order to adapt the previous basis Φlto the Dirichlet boundary conditions Φl(±L)=0,we multiply the previous polynomials by a scale factor. As for the Hermite spectral method(see Section 3.2.1),we use a method of weighted residuals(MWR,see,e.g.,[10,20])and relations(3.17)-(3.20)to obtain the following generalized eigenvalue problem:

whereUNandVNare the vectors containing respectively the coefficientsandl=0,...,N,ofandis the square matrix of order(2N+2):

andBNis the square matrix of order(2N+2):

Hereand=0,...,N,with

As forB0NuN=(uN,Φl′).Note thatis a symmetric matrix with seven diagonal andis a pentadiagonal symmetric matrix.As forB0N=(Φl,Φl′)forlandl′=0,...N(see(3.17)).

To obtain the eigenvalues of the generalized eigenvalue problem(3.21)we use the function DGGEV of the LAPack library which is based on the generalized Schur factorization.

Figure 7: Computation of the eigenvalues for N=50,L=10 and a=1 with the Legendre spectral method.Comparison with the spectral Hermite method and the finite difference method is done.

In Figure 7 we present the solutionsλNof(3.21),computed withN=50,L=10 anda=1.Comparison with the spectral Hermite method(Figure 1)and the finite difference method(Figure 3)is done.We can see that the numerical results are quite similar.

3.2 Another quadratic operator

In this section we consider the following operator:

Fork=2 we retrieve the operator(1.4)studied in the previous section.

We discretize the problemLu=0 using some techniques similar to finite difference methods,with a spatial step equal to one.For simplicity reasons we need to add either periodic boundary conditions or homogeneous Dirichlet boundary conditions.Also we replace Δu(n)byδδ?where:

i.e.

So we have:

3.2.1 Finite difference method with periodic boundary conditions

In this section we are interested to study the problem(3.23)with periodic boundary conditions.So for someN ∈N,we study the following problem:

This gives the following quadratic eigenvalue problem:

where I is theN×Nidentity matrix andA1,A0are given as follows:

with

We start by computing the eigenvalues for different values ofNand for the operatorL.Then,we compute the eigenvalues for some perturbations of the operatorL,i.e.we study the discrete operator:

Figure 8:Eigenvalues of the matrix Ac for N=100,c=1 and k=2.This figure represents a zoom for the case c=1.

for 0≤c≤1 with the same previous periodic boundary conditions.For this we consider the linearization system problem in place of the non-linear problem,so we study the spectrum of the linear systemAcU=λUwith:

wherewithvi=λui,i=1,···,N.A0andA1are given in(3.25)and(3.24)respectively.For the computation of the eigenvalues,we use Matlab.

The results obtained forN=100k=2 andc=1 are presented in Figure 8.The associated domain is[0,N].This figure represents a zoom for the casec=1.We note that the imaginary part of the eigenvaluesλilies between 1.38 and 1.42 in the positive part and between?1.42 and?1.38 in the negative part. Starting from a real partλr=576 all the eigenvalues are aligned on a straight parallel to thex?axiswithλi=1.4141 andλi=?1.4141.The results obtained forN=1000k=2 andc=1 are similar.

In Figure 9 we present the numerical results obtained forN=1000,k=4 and 0≤c ≤1. For the casec=0 we have pure imaginary eigenvalues(since in this case we have just a selfadjoint matrix).The positions of eigenvalues for the casesc=0.2,0.4 confirm the theoretical results.For the casesc=0.6,0.8,1,eigenvalues are localized in a sector delimited by an angle with thex?axissmaller than 2π/6.This is not coherent with the theoretical results.

Figure 9:Eigenvalues of the matrix Ac for N=1000,k=4 and c=0,0.2,0.4,0.6,0.8,1. In the first three figures we can see the cases c=0,0.2,0.4.In the last three figures we can see the cases c=0.6,0.8,1.

3.2.2 Finite difference method with homogeneous Dirichlet boundary conditions

Now we consider the problem(3.23)with homogeneous Dirichlet boundary conditions.So we study the following problem:

We obtain the following nonlinear eigenvalue problem:

where I is theN×Nidentity matrix andA1,A0are given as follows:

Figure 10:Eigenvalues of the matrix Ac for N=2000,k=4,L=10,c=1.

Figure 11:Eigenvalues of the matrix Ac for N=10000,k=4,L=20,c=1.

where

We start by computing the eigenvalues for different values ofNand for the operatorL.Then we compute the eigenvalues for some perturbations of the operatorL,i.e.we consider the discrete operator:

with 0≤c ≤1 and the same previous homogeneous Dirichlet boundary conditions.We do this considering the linearization system problem in place of the non-linear problem.So we study the spectrum of the linear systemAcU=λUwith:

wherewithvi=λui,i=1,···,N.A0andA1are given in(3.28)and(3.27)respectively.We compute the eigenvalues using Matlab.

For the numerical simulations we have considered a domain[?L,+L]and a spatial step Δx=2L/N.For the casek=4,the results obtained forL=10,N=2000(resp.L=20,N=10000)andc=1 are presented in Figures 10 and 11 respectively.For the casek=6,the numerical results obtained for the example(3.26b)withN=10000,c=1 andL=20(resp.L=10)are presented in Figure 12 and 13 respectively.

Figure 12:Eigenvalues of the matrix Ac for N=10000,k=6,L=20,c=1.The figure on the left corresponds to 0≤?λ≤10000 and ?10000≤?λ≤10000.The two figures on the right correspond,up to 0≤?λ≤10000 and 0≤?λ≤10000,down to 0≤?λ≤10000 and ?10000≤?λ≤0.

Figure 13:Eigenvalues of the matrix Ac for N=10000,k=6,L=10,c=1.The figure on the left corresponds to 0≤?λ≤1500 and ?1500≤?λ≤1500.The two figures on the right correspond,up to 0≤?λ≤1500 and 0≤?λ≤1500,down to 0≤?λ≤1500 and ?1500≤?λ≤0.

4 Conclusions and open problems

In this work we have presented a review of some theoretical results obtained for quadratic family of operators:

whereL0andL1are operators in an Hilbert space.

Then we have presented numerical methods to compute the spectrum of such operators.We reduce it to a non self-adjoint linear eigenvalue problem.The numerical methods proposed are based on spectral methods and finite difference methods,in bounded and unbounded domains.For bounded domain we consider homogeneous Dirichlet boundary conditions and periodic boundary conditions.Comparison of the results obtained in unbounded and bounded domains are done.They are based on the size of the containment domain,deduces from the zeros of the Hermite functions.

The numerical results obtained are presented.In particular the numerical instabilities are highlighted. Comparisons of the numerical results obtained in Section 2,with the theoretical results presented in Section 2,are done.These comparisons show the difficulties for the numerical computation of the spectra of such operators.Elimination of the spectral pollution,using staggered grids,and the computation of pseudospectra allow to obtain numerical results in agreement with theoretical results.

A future step in this work is the extension to higher dimension(two and three dimensional case).Indeed,in the multidimensional case very few results are known on the location of the eigenvalues.Numerical simulations,using numerical approaches developed and validate in this article,and completed with parallel computing,should allow to locate eigenvalues.This work is in progress and will be presented elsewhere.

Acknowledgments

This work was initiated during the visit of Fatima Aboud,at the Laboratoire de Mathématiques Jean Leray,Université de Nantes(France),CNRS UMR 6629.This visit was supported by the research projectDéfiMathsof the Fédération de Mathématiques des Pays de la Loire,CNRS FR 2962.F.A.was also supported by the CIMPA(International Center of Pure ans Applied Mathematics).Computations are done thanks to the computer of the CCIPL(Centre de Calcul Intensif des Pays de la Loire).

Appendix A Hermite spectral approximation

The basis{φk}k∈Nof Hermite functions is obtained as an orthonormal basis ofL2(R)of the eigenfunctions of the harmonic oscillator:

We recall briefly its construction(see the basic books of quantum mechanics).Define the creation operatora?and the annihilation operatora:

We satisfy:

Starting by the normalized Gaussian:

verifiedaφ0=0 and thenone define by induction for integerkthe sequence

We verify the following relation by using an algebraic calculation:

where〈,〉denoted the scalar product in the(complex)Hilbert spaceL2(R).We then show thatis a Hilbertian basis ofL2(R).

We have used the following convention:when any integer becomes negative then we replace it by 0.

Bmis equal to the domain ofand the scalar product is equivalent to:

We deduce the following characterization ofBmwith the Hermite coefficient ofu,αk(u):=〈φk,u〉.

Proposition A.1.u∈Bm if and only if

In addition,the scalar product is expressed as the following:

The proposition can be summarized by saying thatBmis identical to the domain of the operatorBy complex interpolation we deduce the intermediate spacesBsfor allspositive reals hence by the duality forsnegative reals.The arguments are identical to the case of usual Sobolev spaces.Fors＜0 theBsare the spaces of temperate distribution.

Hence ifu∈Bmwe have:

More generally we can estimate the error in the spacesBs:

It may be useful to have such Sobolev inequalities explaining the regularity and decay at infinity ofu∈Bsas soon assis large enough:

Proposition A.2.Let m∈N.There exists constants Cm ＞0,Cs,m(m＜2s?2)such that

In particular if m is known and if s ＞m+2then all u in Bs are of class Cm onRand verify the inequality(A.6).

Thus we see that the functionsu ∈Bmare both regular and decreasing to 0 at the infinity more rapidly whenmis big(positive).

Appendix B Pseudospectra

The eigenvalues of Schr?dinger pencils are very unstable(see Section 2).As proposed some times ago by Trefethen[37]it is useful to replace the spectra of non-self adjoint operators by something more stable which is called the pseudospectra.

LetAbe closed operator in the Hilbert spaceHwith domainD(A)dense inH.Recall thatD(A)is an Hilbert space for the graph norm

Definition B.1.The complex number z is in resolvent set ρ(A)of A if and only if A?zIis invertible from D(A)into H and(A?zI)?1∈L(H)where L(H)is the Banach space of linear and continuous maps in H.The spectrum σ(A)

Definition B.2.Consider ε ＞0.The ε-spectrum σε(A)of A is defined as follows.A complex number z∈σε(A)if and only if z∈σ(A)=A?zI.

Proposition B.1.For any matrix A,z ∈σε(A)if and only if smin(A?z)＜ε,where we have denoted smin(A):=min[s(A)].