<figure><img src="image_1.png" /></figure>

HAL
open science

# Numerical approaches to compute spectra of non-self adjoint operators in dimensions two and three

Fatima Aboud, François Jauberteau, Didier Robert

► To cite this version:

Fatima Aboud, François Jauberteau, Didier Robert. Numerical approaches to compute spectra of non-self adjoint operators in dimensions two and three. 2024. hal-04829735

HAL Id: hal-04829735

https://hal.science/hal-04829735v1

Preprint submitted on 10 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Numerical approaches to compute spectra of non-self adjoint operators in dimensions two and three

Fatima Aboud

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

Fatima.Aboud@uodiyala.edu.iq

François Jauberteau & Didier Robert

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

francois.jauberteau@univ-nantes.fr, didier.robert@univ-nantes.fr

## Abstract

In this article we are interested for the numerical computation of spectra of non-self adjoint quadratic operators, in two and three spatial dimensions. Indeed, in the multidimensional case very few results are known on the location of the eigenvalues. This leads to solve nonlinear eigenvalue problems. In introduction we begin with a review of theoretical results and numerical results obtained for the one dimensional case. Then we present the numerical methods developed to compute the spectra (finite difference discretization) for the two and three dimensional cases. The numerical results obtained are presented and analyzed. One difficulty here is that we have to compute eigenvalues of strongly non-self-adjoint operators which are unstable. This work is in continuity of a previous work in one spatial dimension [3].

**Keywords** : non-self adjoint quadratic operators, nonlinear eigenvalue problems, spectra, finite difference methods.

## Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Homogeneous Dirichlet boundary conditions</td><td>3</td></tr><tr><td>2.1</td><td>The two dimensional case</td><td>3</td></tr><tr><td>2.2</td><td>The three dimensional case</td><td>4</td></tr><tr><td>3</td><td>Periodic case</td><td>8</td></tr><tr><td>3.1</td><td>The two dimensional case</td><td>8</td></tr><tr><td>3.2</td><td>The three dimensional case</td><td>9</td></tr><tr><td>4</td><td>Numerical results and discussion</td><td>11</td></tr><tr><td>4.1</td><td>Homogeneous boundary conditions</td><td>11</td></tr><tr><td>4.2</td><td>Periodic boundary conditions</td><td>11</td></tr><tr><td>5</td><td>Conclusion</td><td>13</td></tr><tr><td>6</td><td>Figures</td><td>16</td></tr></table>
# 1 Introduction

The aim of this paper is the extension of [3] (one dimensional in space) to higher dimensional spaces.

In this work we study the polynomial family of operators $L(\lambda) = H_0 + \lambda H_1 + \dots + \lambda^{m-1} H_{m-1} + \lambda^m \mathbb{I}$, where the coefficients $H_0, H_1, \dots, H_{m-1}$ are operators defined on some Hilbert space $\mathcal{H}$, $\mathbb{I}$ is the identity operator and $\lambda$ is a complex parameter. We are interested to study the spectrum of the family $L(\lambda)$.

The problem $L(\lambda)u(x) = 0$, is called a non-linear eigenvalue problem for $m \ge 2$. The number $\lambda_0 \in \mathcal{C}$ is called an eigenvalue of $L(\lambda)$, if there exists $u_0 \in \mathcal{H}$, $u_0 \neq 0$ such that $L(\lambda_0)u_0 = 0$. We consider here a quadratic family ($m = 2$) and in particular we are interested in the case $L_P(\lambda) = \Delta + (P(x)\lambda)^2$, which is defined on the Hilbert space $\mathcal{L}^2(\mathbb{R}^n)$, where $P$ is an elliptic positive polynomial of degree $M \ge 2$. For this example the known results for the existence of eigenvalues are given for $n = 1$ and $n$ is even (see [17]).

The main goal of the work of [1] and [2] is to check the following conjecture, stated by Hellfer-Robert-Wang [17] : for every dimension $n$, for every $M \ge 2$, the spectrum of $L_P$ is non empty, which is confirmed for the following cases :

* $n = 1, 3$, for every polynomial $P$ of degree $M \ge 2$;

* $n = 5$, for every convex polynomial $P$ satisfying some technical conditions;

* $n = 7$, for every convex polynomial $P$.

This result extends to the case of quasi-homogeneous polynomial and quasi-elliptic, for example $P(x, y) = x^2 + y^4$, $x \in \mathbb{R}^{n_1}$, $y \in \mathbb{R}^{n_2}$, $n_1 + n_2 = n$, and $n$ is even. These last results were proved by computing the coefficients of a semi-classical trace formula and by using the theorem of Lidskii.

In [3], some theoretical results are presented for $n = 1$, for quadratic family of operators :

$$L(\lambda) = L_0 + \lambda L_1 + \lambda^2 \mathbb{I}$$

where $L_0$ and $L_1$ are operators on an Hilbert space and $\mathbb{I}$ is the identity operator. Moreover some numerical methods are presented to compute the spectrum of such operators with a non self-adjoint linear eigenvalue problem, in one dimension in space. The proposed numerical methods are based on spectral methods and finite difference methods, in bounded and unbounded domains. For bounded domain homogeneous Dirichlet boundary conditions and periodic boundary conditions are considered with a comparison of the results obtained in unbounded and bounded domains. These comparisons show the difficulties for the numerical computation of the spectra of such operators. The numerical results obtained are in agreement with theoretical results on the location of the eigenvalues (see [10]).

Here we are interested in the location in the complex plane of the eigenvalues for non-linear eigengenvalue problems in two and three dimensional cases. Indeed, in the multidimensional case very few theoretical results are known on the location of the eigenvalues for non-self adjoint quadratic operators. So numerical methods are developed to compute the spectra and numerical results are presented.

Here we notice some references of papers concerning computations of eigenvalues, connected with our paper. For estimation of the spectra of linear differential operators we can refer to [8], [7], [22], [12], [20], [4], [16] and [11]. For numerical methods to solve nonlinear eigenvalue problems (finite dimensional case) and for linearization of such nonlinear eigenvalue problems we can refer to [5], [23], [21], [9], [6], [26], [14] and [15].
# 2 Homogeneous Dirichlet boundary conditions

## 2.1 The two dimensional case

We consider $n = 2$ and the operator

$$L_{\lambda}(u) = -\Delta u + (|x|^2 - \lambda)^2 u$$

where

$$\Delta u(x, y) = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}, (x, y) \in \mathbb{R}^2.$$

We discretize the problem by using finite difference method to solve the following homogeneous problem with the homogeneous Dirichlet boundary conditions

$$\left\{ \begin{array}{l} L_{\lambda}(u(x, y)) = 0, \quad (x, y) \in [-L, L] \times [-L, L] \\ u(L, y) = u(-L, y) = 0 \\ u(x, L) = u(x, -L) = 0 \end{array} \right. \qquad (2.1)$$

So we take the rectangular grid $[-L, L] \times [-L, L]$ with $h = \frac{2L}{N-1}$ for some $N$ ($N$ the number of points on the interval $[-L, L]$). To approximate the Laplacian we need to 5-points $(i, j+1)$, $(i, j-1)$, $(i-1, j)$, $(i+1, j)$ and $(i, j)$ (see Figure 1). For simplicity we denote the point $(x_i, y_j)$ by $(i, j)$ and $u(x_i, y_j)$ by $u_{i,j}$.

We approximate the derivatives $\frac{\partial^2 u}{\partial x^2}$, $\frac{\partial^2 u}{\partial y^2}$ by using Taylor expansion so :

$$\frac{\partial^2 u}{\partial x^2}|_i = \frac{u_{i-1,j} - 2u_{i,j} + u_{i+1,j}}{h^2} + O(h^2)$$

$$\frac{\partial^2 u}{\partial y^2}|_j = \frac{u_{i,j-1} - 2u_{i,j} + u_{i,j+1}}{h^2} + O(h^2)$$

we get

$$\Delta u(x, y)|_{i,j} \approx \frac{u_{i-1,j} - 2u_{i,j} + u_{i+1,j}}{h^2} + \frac{u_{i,j-1} - 2u_{i,j} + u_{i,j+1}}{h^2}$$

$$\Delta u(x, y)|_{i,j} \approx \frac{-4u_{i,j} + u_{i-1,j} + u_{i+1,j} + u_{i,j-1} + u_{i,j+1}}{h^2}.$$

Now we define

$$U_{\bullet,j} = \begin{pmatrix} u_{1,j} \\ u_{2,j} \\ \vdots \\ u_{N,j} \end{pmatrix}_{N \times 1}, \quad \vec{U} = \begin{pmatrix} U_{\bullet,1} \\ U_{\bullet,2} \\ \vdots \\ U_{\bullet,N} \end{pmatrix}_{N^2 \times 1}$$

$$X_{\bullet,j} = \begin{pmatrix} x_1, y_j \\ x_2, y_j \\ \vdots \\ x_N, y_j \end{pmatrix}_{N \times 1}, \quad \vec{X} = \begin{pmatrix} X_{\bullet,1} \\ X_{\bullet,2} \\ \vdots \\ X_{\bullet,N} \end{pmatrix}_{N^2 \times 1}$$

$$|X_{\bullet,j}|^2 = \begin{pmatrix} x_1^2 + y_j^2 \\ x_2^2 + y_j^2 \\ \vdots \\ x_N^2 + y_j^2 \end{pmatrix}_{N \times 1}, \quad |\vec{X}|^2 = \begin{pmatrix} |X_{\bullet,1}|^2 \\ |X_{\bullet,2}|^2 \\ \vdots \\ |X_{\bullet,N}|^2 \end{pmatrix}_{N^2 \times 1}$$
To construct the Laplacian matrix we need to define the matrix $D_4$ which is an $N \times N$ matrix such that

$$D_4 = \begin{pmatrix} -4 & 1 & 0 & \cdots & \cdots & 0 \\ 1 & -4 & 1 & 0 & \cdots & 0 \\ 0 & 1 & -4 & 1 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & \cdots & 0 & 1 & -4 & 1 \\ 0 & \cdots & \cdots & 0 & 1 & -4 \\ 0 & \cdots & \cdots & \cdots & 0 & 1 & -4 \end{pmatrix}_{N \times N}$$

We denote the Laplacian matrix by $\Delta$ which has the following form

$$h^2 \Delta = \begin{pmatrix} D_4 & I & 0 & \cdots & \cdots & 0 \\ I & D_4 & I & 0 & \cdots & 0 \\ 0 & I & D_4 & I & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & \cdots & 0 & I & D_4 & I \\ 0 & \cdots & \cdots & 0 & I & D_4 \\ 0 & \cdots & \cdots & \cdots & 0 & D_4 \end{pmatrix}_{N^2 \times N^2} \quad (2.2)$$

where $I$ is the identity $N \times N$ matrix and we note that the zeros represent the zero $N \times N$ matrices. So we have the following discrete problem

$$- \Delta \vec{U} + (|\vec{X}|^2 - \lambda \mathbb{I})^2 \vec{U} = 0 \quad (2.3)$$

i.e.

$$(-\Delta + \mathrm{diag}(|\vec{X}|^4))\vec{U} - 2\lambda \mathrm{diag}(|\vec{X}|^2) \vec{U} + \lambda^2 \mathbb{I} \vec{U} = 0$$

We set

$$H_0 = -\Delta + \mathrm{diag}(|\vec{X}|^4), \quad H_1 = -2 \mathrm{diag}(|\vec{X}|^2)$$

where $\mathrm{diag}(|\vec{X}|^4)$, $\mathrm{diag}(|\vec{X}|^2)$ represent square diagonal matrices with the elements of vector $|\vec{X}|^4$, $|\vec{X}|^2$ on the main diagonal, respectively. By using the linearization of this problem we get

$$\mathcal{A} = \begin{pmatrix} 0 & I \\ -H_0 & -H_1 \end{pmatrix} \quad (2.4)$$

## 2.2 The three dimensional case

We consider the case of the space dimension $n = 3$ for the family of operators:

$$L_{\lambda} = -\Delta + (|\mathbf{x}|^2 - \lambda)^2 \mathbb{I} = -\Delta + |\mathbf{x}|^4 \mathbb{I} - 2\lambda |\mathbf{x}|^2 \mathbb{I} + \lambda^2 \mathbb{I}$$

with $\mathbf{x} = (x, y, z) \in \mathbb{R}^3$ and $\mathbb{I}$ the identity operator. Introducing a real parameter $c$ we study a more general case of this family of operators

$$L_{\lambda,c} = -\Delta + |\mathbf{x}|^4 - 2c\lambda|\mathbf{x}|^2 + \lambda^2$$

$$\left\{ \begin{array}{l} L_{\lambda,c}(u(x, y, z)) = 0, \quad (x, y, z) \in [-L, L] \times [-L, L] \times [-L, L] \\ u(L, y, z) = u(-L, y, z) = 0 \\ u(x, L, z) = u(x, -L, z) = 0 \\ u(x, y, L) = u(x, y, -L) = 0 \end{array} \right. \quad (2.5)$$
<figure><img src="image_2.png" /><figcaption>Figure 1: 5-points finite difference 2D</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: 7-points finite difference 3D</figcaption></figure>

For simplicity reason, we denote the point $(x_i, y_j, z_k)$ by $(i, j, k)$ and $u(x_i, y_j, z_k)$ by $u_{i,j,k}$.
To calculate the Laplacian we need to 7-points $(i, j, k)$, $(i-1, j, k)$, $(i+1, j, k)$, $(i, j-1, k)$, $(i, j+1, k)$, $(i, j, k-1)$, $(i, j, k+1)$ (see Figure 2). We have :

$$ \Delta u(x, y, z) = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2}, \quad (x, y, z) \in \mathbb{R}^3 $$

and we approximate the derivatives $\frac{\partial^2 u}{\partial x^2}$, $\frac{\partial^2 u}{\partial y^2}$ and $\frac{\partial^2 u}{\partial z^2}$ by using Taylor expansion :

$$ \frac{\partial^2 u}{\partial x^2}|_i = \frac{u_{i-1,j,k} - 2u_{i,j,k} + u_{i+1,j,k}}{h^2} + O(h^2) $$

$$ \frac{\partial^2 u}{\partial y^2}|_j = \frac{u_{i,j-1,k} - 2u_{i,j,k} + u_{i,j+1,k}}{h^2} + O(h^2) $$

$$ \frac{\partial^2 u}{\partial z^2}|_k = \frac{u_{i,j,k-1} - 2u_{i,j,k} + u_{i,j,k+1}}{h^2} + O(h^2) $$

So we obtain :

$$ \Delta u(x, y, z)|_{i,j,k} \approx \frac{u_{i-1,j,k} - 2u_{i,j,k} + u_{i+1,j,k}}{h^2} + \frac{u_{i,j-1,k} - 2u_{i,j,k} + u_{i,j+1,k}}{h^2} + \frac{u_{i,j,k-1} - 2u_{i,j,k} + u_{i,j,k+1}}{h^2} $$

i.e.

$$ \Delta u(x, y, z)|_{i,j,k} \approx \frac{-6u_{i,j,k} + u_{i-1,j,k} + u_{i+1,j,k} + u_{i,j-1,k} + u_{i,j+1,k} + u_{i,j,k-1} + u_{i,j,k+1}}{h^2} $$

Now we define

$$ U_{\bullet,j,k} = \begin{pmatrix} u_{1,j,k} \\ u_{2,j,k} \\ \vdots \\ u_{N,j,k} \end{pmatrix}_{N \times 1}, \quad V_{i,\bullet,k} = \begin{pmatrix} U_{\bullet,1,k} \\ U_{\bullet,2,k} \\ \vdots \\ U_{\bullet,N,k} \end{pmatrix}_{N^2 \times 1}, \quad \vec{U} = \begin{pmatrix} V_{i,\bullet,1} \\ V_{i,\bullet,2} \\ \vdots \\ V_{i,\bullet,N} \end{pmatrix}_{N^3 \times 1} $$
$$X_{\bullet,j,k} = \begin{pmatrix} x_1, y_j, z_k \\ x_2, y_j, z_k \\ \vdots \\ x_N, y_j, z_k \end{pmatrix}_{N \times 1}, \quad X'_{i,\bullet,k} = \begin{pmatrix} X_{\bullet,1,k} \\ X_{\bullet,2,k} \\ \vdots \\ X_{\bullet,3,k} \end{pmatrix}_{N^2 \times 1}, \quad \vec{X} = \begin{pmatrix} X'_{i,\bullet,1} \\ X'_{i,\bullet,2} \\ \vdots \\ X'_{i,\bullet,N} \end{pmatrix}_{N^3 \times 1}$$

$$|\vec{X}|^2 = \begin{pmatrix} x_i^2 + y_\bullet^2 + z_1^2 \\ x_i^2 + y_\bullet^2 + z_2^2 \\ \vdots \\ x_i^2 + y_\bullet^2 + z_N^2 \end{pmatrix}_{N^3 \times 1}, \quad |\vec{X}|^4 = \begin{pmatrix} (x_i^2 + y_\bullet^2 + z_1^2)^2 \\ (x_i^2 + y_\bullet^2 + z_2^2)^2 \\ \vdots \\ (x_i^2 + y_\bullet^2 + z_N^2)^2 \end{pmatrix}_{N^3 \times 1}$$

To construct the matrix of Laplacian $\Delta$ we need to define the following matrices. Let $D_6$ be an $N \times N$ matrix defined as follows

$$D_6 = \begin{pmatrix} -6 & 1 & 0 & \cdots & \cdots & 0 \\ 1 & -6 & 1 & 0 & \cdots & 0 \\ 0 & 1 & -6 & 1 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & \cdots & 0 & 1 & -6 & 1 \\ 0 & \cdots & \cdots & 0 & 1 & -6 \\ 0 & \cdots & \cdots & \cdots & 0 & 1 & -6 \end{pmatrix}_{N \times N}$$

Now we define $G_6$ which is an $N^2 \times N^2$ matrix such that

$$G_6 = \begin{pmatrix} D_6 & I & 0 & \cdots & \cdots & 0 \\ I & D_6 & I & 0 & \cdots & 0 \\ 0 & I & D_6 & I & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & \cdots & 0 & I & D_6 & I & 0 \\ 0 & \cdots & \cdots & 0 & I & D_6 & I \\ 0 & \cdots & \cdots & \cdots & 0 & I & D_6 \end{pmatrix}_{N^2 \times N^2}$$

where $I$ is the $N \times N$ identity and the zeros represent the $N \times N$ zero matrix. So the 3D Laplacian finite difference matrix is given by

$$h^2 \Delta = \begin{pmatrix} G_6 & J & 0 & \cdots & \cdots & 0 \\ J & G_6 & J & 0 & \cdots & 0 \\ 0 & J & G_6 & J & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & \cdots & 0 & J & G_6 & J & 0 \\ 0 & \cdots & \cdots & 0 & J & G_6 & J \\ 0 & \cdots & \cdots & \cdots & 0 & J & G_6 \end{pmatrix}_{N^3 \times N^3} \quad (2.6)$$

where $J$ is the $N^2 \times N^2$ identity and the zeros are the zero $N^2 \times N^2$ matrix. We note that the Laplacian matrix will be $N^3 \times N^3$.

So we have the following discrete problem

$$- \Delta \vec{U} + (|\vec{X}|^2 - \lambda I)^2 \vec{U} = 0 \quad (2.7)$$
With the parameter $c$ we deduce the more general discrete problem

$$(-\Delta + \mathrm{diag}(|\vec{X}|^4))\vec{U} - 2\lambda c \mathrm{diag}(|\vec{X}|^2) \vec{U} + \lambda^2 \mathbb{I} \vec{U} = 0$$

We set

$$H_0 = -\Delta + \mathrm{diag}(|\vec{X}|^4), \quad H_1 = -2 c \mathrm{diag}(|\vec{X}|^2)$$

By using the linearization of this problem we get

$$A = \begin{pmatrix} 0 & \mathbb{I} \\ -H_0 & -H_1 \end{pmatrix} \qquad (2.8)$$

with

$$H_0 = \begin{pmatrix} \Theta_{i, \bullet, 1} - \mathbb{G}_6/h^2 & -\mathbb{J} & 0 & \dots & \dots & 0 \\ -\mathbb{J} & \Theta_{i, \bullet, 2} - \mathbb{G}_6/h^2 & -\mathbb{J} & 0 & \dots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & \dots & 0 & -\mathbb{J} & \Theta_{i, \bullet, N-1} - \mathbb{G}_6/h^2 & -\mathbb{J} \\ 0 & \dots & \dots & 0 & -\mathbb{J} & \Theta_{i, \bullet, N} - \mathbb{G}_6/h^2 \end{pmatrix}_{N^3 \times N^3}$$

where $\Theta_{i, \bullet, j} = \mathrm{diag}(|X_{i, \bullet, j}|^4)$ is the diagonal $N^2 \times N^2$ matrix with $|X'_{i, \bullet, j}|^4$ on the diagonal, and $H_1 = -2c \mathrm{diag}(|\vec{X}|^2)$ with $\mathrm{diag}(|\vec{X}|^2)$ is the diagonal $N^3 \times N^3$ matrix with $|\vec{X}|^2$ on the diagonal.
## 3 Periodic case

Consider the following differential operator :

$$ - \Delta u + \lambda^2 u + \frac{\lambda}{i} \sum_{j=1}^{n} a_j(x_j) \frac{\partial u}{\partial x_j} = 0, \quad x \in \mathbb{R}^n \qquad (3.9) $$

with periodic boundary conditions.

In the case where the coefficients $a_j(x_j)$ are space independent, if we consider $u = \sum_{\mathbf{k} \in \mathbb{Z}} \hat{u}_{\mathbf{k}} e^{i \mathbf{k} \cdot \mathbf{x}}$ and substituting in (3.9) we obtain :

$$ k^2 \hat{u}_{\mathbf{k}} + \lambda^2 \hat{u}_{\mathbf{k}} + \frac{\lambda}{i} \sum_{j=1}^{n} ia_j k_j \hat{u}_{\mathbf{k}} = 0 \\ \Rightarrow \lambda^2 + \lambda \sum_{j=1}^{n} a_j k_j + k^2 = 0 $$

where $\mathbf{x} = (x_i)_{i=1,...,n}$, $\mathbf{k} = (k_i)_{i=1,...,n}$ and $k = |\mathbf{k}|$.

In the following we consider the case where the coefficients are space dependent and periodic in space.

### 3.1 The two dimensional case

Let $\Omega = (0, L) \times (0, L) \subset \mathbb{R}^2$, we study the following problem

$$ -\Delta u + \lambda^2 u + \frac{\lambda}{i} \left( a_1(x) \frac{\partial u}{\partial x} + a_2(y) \frac{\partial u}{\partial y} \right) = 0 \\ u(x+L, y+L) = u(x, y) $$

where $L \in \mathbb{R}$, $x \in [0, L]$ and $y \in [0, L]$. So we take the rectangular grid $[0, L] \times [0, L]$ with $h = \frac{2L}{N-1}$ for some $N$ ($N$ the number of points on the interval $[0, L]$). We set :

$$ H_0 = -\Delta, \quad H_1 = \frac{1}{i} \sum_{j=1}^{2} a_j(x_j) \frac{\partial}{\partial x_j} $$

To get the formula of the linearization operator $\mathcal{A}$ we need to have the formulas of the matrices $H_0, H_1$. We note that $H_0$ is the Laplacian, so it has the formula given in (2.2). To calculate the Jaccobi matrix we need to 4-points $(i, j+1), (i, j-1), (i-1, j), (i+1, j)$ (for simplicity we denote the point $(x_i, y_j)$ by $(i, j)$ and $u(x_i, y_j)$ by $u_{i,j}$). We approximate the derivatives $\frac{\partial u}{\partial x}, \frac{\partial u}{\partial y}$ by using Taylor expansion so

$$ \frac{\partial u}{\partial x}|_i = \frac{u_{i+1,j} - u_{i-1,j}}{2h} + O(h^2) $$

$$ \frac{\partial u}{\partial y}|_j = \frac{u_{i,j+1} - u_{i,j-1}}{2h} + O(h^2) $$
To give the formula of $H_1$ we need the following matrices

$$ \mathbb{D}_x = \begin{pmatrix} 0 & a_1(x_1) & 0 & \cdots & \cdots & -a_1(x_1) \\ -a_1(x_2) & 0 & a_1(x_2) & 0 & \cdots & 0 \\ 0 & -a_1(x_3) & 0 & a_1(x_3) & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & \cdots & 0 & -a_1(x_{N-2}) & 0 & a_1(x_{N-2}) & 0 \\ 0 & \cdots & \cdots & 0 & -a_1(x_{N-1}) & 0 & a_1(x_{N-1}) \\ a_1(x_N) & \cdots & \cdots & \cdots & 0 & -a_1(x_N) & 0 \end{pmatrix}_{N \times N} $$

so

$$ H_1 = \frac{1}{2ih} \begin{pmatrix} \mathbb{D}_x & Dy_1 & 0 & \cdots & \cdots & -Dy_1 \\ -Dy_2 & \mathbb{D}_x & Dy_2 & 0 & \cdots & 0 \\ 0 & -Dy_3 & \mathbb{D}_x & Dy_3 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & \cdots & 0 & -Dy_{N-2} & \mathbb{D}_x & Dy_{N-2} & 0 \\ 0 & \cdots & \cdots & 0 & -Dy_{N-1} & \mathbb{D}_x & Dy_{N-1} \\ Dy_N & \cdots & \cdots & \cdots & 0 & -Dy_N & \mathbb{D}_x \end{pmatrix}_{N^2 \times N^2} \quad (3.10) $$

where $Dy_j = \text{diag}(a_2(y_j))$.

## 3.2 The three dimensional case

Let $\Omega = (0, L) \times (0, L) \times (0, L) \subset \mathbb{R}^3$, we study the following problem

$$ -\Delta u + \lambda^2 u + \frac{\lambda}{i} \left( a_1(x) \frac{\partial u}{\partial x} + a_2(y) \frac{\partial u}{\partial y} + a_3(z) \frac{\partial u}{\partial z} \right) = 0 \\ u(x+L, y+L, z+L) = u(x, y, z) $$

where $x \in [0, L]$, $y \in [0, L]$, $z \in [0, L]$ and $L \in \mathbb{R}$. So we take the rectangular grid $[0, L] \times [0, L] \times [0, L]$ with $h = \frac{2L}{N-1}$ for some $N$ ($N$ the number of points on the interval $[0, L]$). We set :

$$ H_0 = -\Delta, \quad H_1 = \frac{1}{i} \sum_{j=1}^{3} a_j(x_j) \frac{\partial}{\partial x_j} $$

To get the formula of the linearization operator $\mathcal{A}$ we need to have the formulas of the matrices $H_0, H_1$. We note that $H_0$ is the Laplacian, so it has the formula given in (2.6).

For simplicity reason, we denote the point $(x_i, y_j, z_k)$ by $(i, j, k)$ and $u(x_i, y_j, z_k)$ by $u_{i,j,k}$.

To calculate the Jaccobi matrix we need to 6-points $(i-1, j, k), (i+1, j, k), (i, j-1, k), (i, j+1, k), (i, j, k-1), (i, j, k+1)$. We approximate the derivatives $\frac{\partial u}{\partial x}$, $\frac{\partial u}{\partial y}$ and $\frac{\partial u}{\partial z}$ by using Taylor expansion so :

$$ \frac{\partial u}{\partial x}|_i = \frac{u_{i+1,j,k} - u_{i-1,j,k}}{2h} + O(h^2) $$

$$ \frac{\partial u}{\partial y}|_j = \frac{u_{i,j+1,k} - u_{i,j-1,k}}{2h} + O(h^2) $$

$$ \frac{\partial u}{\partial z}|_k = \frac{u_{i,j,k+1} - u_{i,j,k-1}}{2h} + O(h^2) $$
To give the formula of $H_1$, we need the following matrices

$$ \mathbb{D}_x = \begin{pmatrix} 0 & a_1(x_1) & 0 & \cdots & \cdots & -a_1(x_1) \\ -a_1(x_2) & 0 & a_1(x_2) & 0 & \cdots & 0 \\ 0 & -a_1(x_3) & 0 & a_1(x_3) & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & \cdots & 0 & -a_1(x_{N-2}) & 0 & a_1(x_{N-2}) & 0 \\ 0 & \cdots & \cdots & 0 & -a_1(x_{N-1}) & 0 & a_1(x_{N-1}) \\ a_1(x_N) & \cdots & \cdots & \cdots & 0 & -a_1(x_N) & 0 \end{pmatrix}_{N \times N} $$

$$ Dy_j = (diag(a_2(y_j)))_{N \times N} \quad \text{and} \quad Dz_j = (diag(a_3(z_j)))_{N^2 \times N^2} $$

so

$$ \mathbb{G}_{x,y} = \begin{pmatrix} \mathbb{D}_x & Dy_1 & 0 & \cdots & \cdots & -Dy_1 \\ -Dy_2 & \mathbb{D}_x & Dy_2 & 0 & \cdots & 0 \\ 0 & -Dy_3 & \mathbb{D}_x & Dy_3 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & \cdots & 0 & -Dy_{N-2} & \mathbb{D}_x & Dy_{N-2} & 0 \\ 0 & \cdots & \cdots & 0 & -Dy_{N-1} & \mathbb{D}_x & Dy_{N-1} \\ Dy_N & \cdots & \cdots & \cdots & 0 & -Dy_N & \mathbb{D}_x \end{pmatrix}_{N^2 \times N^2} \quad (3.11) $$

thus

$$ H_1 = \frac{1}{2ih} \begin{pmatrix} \mathbb{G}_{x,y} & Dz_1 & 0 & \cdots & \cdots & -Dz_1 \\ -Dz_2 & \mathbb{G}_{x,y} & Dz_2 & 0 & \cdots & 0 \\ 0 & -Dz_3 & \mathbb{G}_{x,y} & Dz_3 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & \cdots & 0 & -Dz_{N-2} & \mathbb{G}_{x,y} & Dz_{N-2} & 0 \\ 0 & \cdots & \cdots & 0 & -Dz_{N-1} & \mathbb{G}_{x,y} & Dz_{N-1} \\ Dz_N & \cdots & \cdots & \cdots & 0 & -Dz_N & \mathbb{G}_{x,y} \end{pmatrix}_{N^3 \times N^3} \quad (3.12) $$
# 4 Numerical results and discussion

In this section we present the numerical results obtained in dimensions 2 and 3 for homogeneous boundary conditions and for periodic boundary conditions. The numerical computation of the spectra of non-self adjoint operators being unstable, it is possible to compute the pseudospectra, see for example [13], [18], [19], [24], [25]. Indeed the computation of the pseudospectra is more stable, but it is expensive since it requires to compute the minimal singular value $s_{min}(A - zI)$ at each point $z$ of the mesh. Here we have computed the spectra using Lapack library.

## 4.1 Homogeneous boundary conditions

In this section we present the numerical results obtained for homogeneous boundary conditions and spatial dimensions $n = 2$ and $n = 3$. The operator considered here is $L_{\lambda,c}$ defined in (2.5).

We start by presenting the results for $c = 0$ and $N = 100$ (resp. $N = 25$) for the two (resp. three) dimensional case. In this case, the spectrum contains just a purely imaginary eigenvalues, see Figure 3.

When the parameter $c$ is increased ($c = 0.5, 1, 2, 5$ and $10$) the real part of the eigenvalues $\lambda$ is increased, see Figures 4–8. In one spatial dimension and $c = 1$ according the theory we have that $Sp(L_{\lambda,c})$ is included in the two sectors $\{\lambda \in \mathbb{C}, |\arg(\lambda)| \ge \frac{\pi}{3}\}$ (see [10] and [3] for a review of theoretical results). Numerical results obtained here show that this location of the eigenvalues is still satisfied in dimensions 2 and 3.

The CPU time is expensive in 2D and in 3D, since the size of the matrices obtained after linearization and discretization is $2N^n \times 2N^n$, with $n = 2$ and $n = 3$. So the spatial resolution is quite limited. We have studied the influence of the spatial resolution on the numerical results for $c = 1$. On Figures 9a and 10a we test the influence of the choice of the parameter $N$ with a constant size of the domain $[-L, L]$. We can observe a similar comportment in 2D and in 3D. When $N$ is increased the repartition of the imaginary part of the eigenvalues is increased but it is not the case for the real part. However on Figures 9b and 10b when the resolution $N$ is constant and the size of the domain is increased ($L = 0.5, 1$ and $2$) we can see that the repartition of the real and imaginary parts of the eigenvalues are increased. So, in conclusion increasing the parameter $N$ has an effect on the repartition of the imaginary part. Moreover, increasing the size of the domain has an effect on the repartition on the real and imaginary parts of the eigenvalues.

In conclusion, increasing $c$ increased the repartition of the spectrum along the real part. Increasing $N$ increased the repartition of the spectrum along the imaginary part. And increasing $L$ increased the repartition of the spectrum along the real and imaginary parts.

## 4.2 Periodic boundary conditions

In this section we present the numerical results obtained for periodic boundary conditions and spatial dimensions $n = 2$ and $n = 3$. The operator considered here is defined in (3.9). The domain retained for the simulations is $[-\pi, \pi]$. For the spatial resolution parameter $N$ we have retained $N = 100$ (resp. $N = 25$) for the simulations in one (resp. two) dimension in space. And we consider different choices for the coefficients $a_j(x_j)$, $j = 1, \dots, n$, periodic in space.

Firstly we study the two dimensional case with constant coefficients $a_j(x_j)$, $j = 1, 2$. We have presented the results obtained with different choices of the coefficients : $a_1(x_1) = a_2(x_2) = 1$ on Figure 11a, $a_1(x_1) = 1$, $a_2(x_2) = \sqrt{2}$ on Figure 11b and $a_1(x_1) = 1$, $a_2(x_2) = 5\sqrt{2}$ on Figure
12a. We can observe that the repartition of the real (resp. imaginary) parts of the eigenvalues are symetric with respect to the vertical (resp. horizontal) axes, which was not the case for the previous operator (2.5). Moreover if the repartition of the eigenvalues is quite similar on Figures 11a and 11b (the real part increases with the imaginary part in modulus), it is different for the eigenvalues with small imaginary part on Figure 12a (the real part decreases when the imaginary part increases in modulus).

Then we consider the two dimensional case with sinusoidal coefficients $a_1(x_1) = \sin(x_1)$ and $a_2(x_2) = \sin(x_2)$. The numerical results are presented on Figure 12b. We can see that the repartition of the eigenvalues is always symetric, as for constant coefficients. However the real part of the eigenvalues are smaller in modulus than for constant coefficients. Moreover the real part of the eigenvalues is approximatively equal to zero when the modulus of the imaginary part is less than 20.

Now we present the three dimensional case with different choices of the coefficients. Numerical results for constant coefficients $a_j(x_j) = 1, i = 1, \dots, 3$ and $a_1(x_1) = 1, a_2(x_2) = \sqrt{2}, a_3(x_3) = 5\sqrt{2}$ are represented on Figures 13a and 13b respectively. Then we consider sinusoidal coefficients $a_1(x_1) = \sin(x_1)$, $a_2(x_2) = \sin(x_2)$ and $a_3(x_3) = \sin(x_3)$. The numerical results are presented on Figure 14a. As for the periodic two dimensional case, we can see that the repartition of the eigenvalues is symetric with respect the horizontal and vertical axes when the coefficients $a_j(x_j) = 1, i = 1, \dots, 3$ are the same (see Figures 13a and 14a). However, when these coefficients are different these symetries are lost for the eigenvalues with small imaginary part (symetry is with respect the origin, see Figure 13b). Moreover, apart from the eigenvalues with small imaginary parts, the eigenvalues are localized on the imaginary axis.
# 5 Conclusion

In this work we have considered quadratic family of operators :

$$L(\lambda) = L_0 + \lambda L_1 + \lambda^2 \mathbb{I}$$

where $L_0$ and $L_1$ are operators on an Hilbert space and $\mathbb{I}$ is the identity operator.

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 finite difference methods, in bounded domains. We consider homogeneous Dirichlet boundary conditions and periodic boundary conditions.

The numerical results obtained in the two and three dimensional cases are presented since very few theoretical results are known on the location of the eigenvalues in the multidimensional case. Numerical results obtained here show that theoretical results on the location of the eigenvalues in one dimension in space are still satisfied in dimensions 2 and 3.

# Acknowledgments

This work has benefit of the visit of Fatima Aboud, at the Laboratoire de Mathématiques Jean Leray, Nantes Université (France). This visit was supported by the Fédération de Mathématiques des Pays de la Loire, CNRS FR 2962.
References

[1] F. Aboud. Ph.D Thesis, *Problèmes aux valeurs propres non-linéaires*. University of Nantes, May 2009. https://tel.archives-ouvertes.fr/tel-00410455.

[2] F. Aboud and D. Robert. Asymptotic expansion for nonlinear eigenvalue problems. Journal de Mathématiques Pures et Appliquées **93** (2), pp.149-162 (2010).

[3] F. Aboud, F. Jauberteau, G. Moebs, D. Robert. Numerical Approaches to Compute Spectra of Non-Self Adjoint Operators and Quadratic Pencils. Journal of Mathematical Study **53** (1), pp.12-44 (2020).

[4] R. Alam. On Spectral Approximation of Linear Operators. Journal of Mathematical Analysis and Applications **226**, pp.229-244 (1998).

[5] M. Asch and G. Lebeau. The spectrum of the damped wave operator for a bounded domain in $\mathbb{R}^2$, Experimental Mathematics. **12** (2), pp.227-241 (2003).

[6] D. Bindel and A. Hood. Localization theorems for nonlinear eigenvalue problems. SIAM Journal on Matrix Analysis and Applications **34** (4), pp.1728-1749 (2013).

[7] J.H. Bramble and J.E. Osborn. Rate of convergence estimates for nonselfadjoint eigenvalue approximations. Mathematics of Computation **27**, pp.525-549 (1973).

[8] J.H. Bramble and A.H. Schatz. Rayleigh-Ritz-Galerkin methods for Dirichlet's problem using subspaces without boundary conditions. Comm. Pure Appl. Math. **23**, pp.653-675 (1970).

[9] E. Cancès, R. Chakir and Y. Maday. Numerical analysis of nonlinear eigenvalue problems. J. Sci. Comput. **45** (2), pp.90-117 (2010).

[10] M. Christ. Analytic hypoellipticity, representation of nilpotent groups and non-linear eigenvalue problem. Duke Math Journal **72** (3), pp.595-639 (1993).

[11] F. Chatelin. *Spectral approximation of linear operators*. Classics in applied mathematics **65** Society for Industrial and Applied Mathematics, Philadelphia (2011).

[12] F. Chatelin and J. Lemordant. Error bounds in the approximation of eigenvalues of differential and integral operators. Journal of Mathematical Analysis and Applications **62**, pp.257-271 (1978).

[13] E.B. Davies. Pseudospectra, the harmonic oscillator and complex resonances. R. Soc. Lond. Proc. Ser. A Math. Phys. EngSci. **455**, pp.585-599 (1999)

[14] C. Effenberger. Ph.D Thesis, *Robust solution methods for nonlinear eigenvalue problems*. École Polytechnique Fédérale de Lausanne, Août 2013.

[15] L. Grammont, N.J. Higham and F. Tisseur. A framework for analyzing nonlinear eigenproblems and parametrized linear systems. Linear Algebra and its Applications **435**, pp.623-640 (2011).

[16] A.C. Hansen. On the approximation of spectra of linear operators on Hilbert spaces. Journal of functional analysis **254**, pp.2092-2126 (2008).

[17] B. Helffer, D. Robert and Xue Ping Wang. Semiclassical analysis of a nonlinear eigenvalue problem and analytic hypoellipticity. Algebra i Analiz **16** (1), pp. 285-296 (2004).

[18] N. Higham and F. Tisseur. Structured pseudo spectra for polynomial eigenvalue problems, with applications. SIAM. Matrix Anal. Appl. **23** (1), pp.187-208 (2001).
[19] N. Higham and F. Tisseur. More on pseudo spectra for polynomial eigenvalue problems and application in control theory. Linear Algebra and its applications **351-352**, pp.435-453 (2002).

[20] W. Kolata and J. Osborn. Non-selfadjoint spectral approximation and the finite element method. Springer Verlag, Lecture Notes in Mathematics **701**, pp.115-133 (1979).

[21] V. Mehrmann and H. Voss. Nonlinear eigenvalue problems : a challenge for modern eigenvalue methods, Mitt. der Ges. für Angwandte Mathematik and Mechanik **27**, pp.121-151 (2005).

[22] J.E. Osborn. Spectral approximation for compact operators. Mathematics of Computation **29** (131), pp.712-725 (1975).

[23] F. Tisseur and K. Meerbergen. The quadratic eigenvalue problem. SIAM Review **43**, pp.235-286 (2001).

[24] L.N. Trefethen. Computation of pseudospectra. Acta Numer. **8**, pp.247-295 (1999).

[25] L.N. Trefethen and M. Embree. *Spectra an pseudospectra*. Princeton University Press (2005).

[26] R. Van Beeumen, K. Meerbergen and W. Michiels. A rational Krylov method based on Hermite interpolation for nonlinear eigenvalue problems. SIAM J. Sci. Comput. **35** (1), pp.A327-A350 (2013).
# 6 Figures

<figure><img src="image_4.png" /><figcaption>Figure 3: Homogeneous boundary conditions for $n = 2, 3$ and $c = 0$</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Homogeneous boundary conditions for $n = 2, 3$ and $c = 0.5$</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Homogeneous boundary conditions for $n = 2, 3$ and $c = 1$</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Homogeneous boundary conditions for $n = 2, 3$ and $c = 2$</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Homogeneous boundary conditions for $n = 2, 3$ and $c = 5$</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8: Homogeneous boundary conditions for $n = 2, 3$ and $c = 10$</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9: Homogenous boundary conditions for $n = 2$ and $c = 1$</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10: Homogenous boundary conditions for $n = 3$ and $c = 1$</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 11: Periodic boundary conditions for $n = 2$</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12: Periodic boundary conditions for $n = 2$</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 13: Periodic boundary conditions for $n = 3$</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 14: Periodic boundary conditions for $n = 3$</figcaption></figure>