The spectral-Galerkin approximation of nonlinear eigenvalue problems

Accepted Manuscript The spectral-Galerkin approximation of nonlinear eigenvalue problems

Jing An, Jie Shen, Zhimin Zhang

PII: DOI: Reference:

S0168-9274(18)30101-6 https://doi.org/10.1016/j.apnum.2018.04.012 APNUM 3366

To appear in:

Applied Numerical Mathematics

Received date: Revised date: Accepted date:

3 June 2017 20 February 2018 19 April 2018

Please cite this article in press as: J. An et al., The spectral-Galerkin approximation of nonlinear eigenvalue problems, Appl. Numer. Math. (2018), https://doi.org/10.1016/j.apnum.2018.04.012

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

The Spectral-Galerkin Approximation of Nonlinear Eigenvalue Problems Jing An

∗

Jie Shen

†

Zhimin Zhang‡

Abstract In this paper we present and analyze a polynomial spectral-Galerkin method for nonlinear elliptic eigenvalue problems of the form −div(A∇u)+V u+f (u2 )u = λu, uL2 = 1. We estimate errors of numerical eigenvalues and eigenfunctions. Spectral accuracy is proved under rectangular meshes and certain conditions of f . In addition, we establish optimal error estimation of eigenvalues in some hypothetical conditions. Then we propose a simple iteration scheme to solve the underlying an eigenvalue problem. Finally, we provide some numerical experiments to show the validity of the algorithm and the correctness of the theoretical results. Keywords: Spectral-Galerkin approximation, error estimation, iteration algorithm, nonlinear eigenvalue problems

1

Introduction

Eigenvalue problems appear in many mathematical models for scientific and engineering applications, such as the calculation of the vibration modes of a mechanical structure in the framework of nonlinear elasticity, the Gross-Pitaevskii equation describing the steady states of Bose-Einstein condensates [14, 2], and the Hartree-Fock and Kohn-Sham equations used to calculate ground state electronic structures of molecular systems in quantum chemistry and materials science [6, 12, 16]. However, most of the existing analysis for eigenvalue approximations are concerned with linear eigenvalue problems [1], and there are relatively few results concerning approximation of nonlinear eigenvalue problems [18, 17, 7, 9, 10, 11, 13], and most of them are based on finite element methods with an exception in [7, 8] where an error estimate for Fourier spectral method to a periodic nonlinear eigenvalue problem is derived. To the best of our knowledge, there has no report on high order numerical methods for non-periodic nonlinear eigenvalue problems. Thus, the aim of this paper is to develop and analyze a spectral Galerkin method for a nonlinear elliptic eigenvalue problem. In particular, we shall extend the error estimates established in [7] for eigenvalues and eigenfunctions for periodic case ∗ Beijing Computational Science Research Center, Beijing 100193, China ([email protected]); and School of Mathematical Sciences, Guizhou Normal University, Guiyang 550025, China ([email protected]). † Fujian Provincial Key Laboratory on Mathematical Modeling and High Performance Scientific Computing and School of Mathematical Science, Xiamen University, Xiamen 361005, PR China; and Department of Mathematics Purdue University, West Lafayette, IN 47907, USA ([email protected]). ‡ Beijing Computational Science Research Center, Beijing 100193, China ([email protected]); and Department of Mathematics, Wayne State University, Detroit, MI 48202, USA ([email protected]).

1

with Fourier approximation to non-periodic case with polynomial approximation. We also describe an efficient implementation of the spectral-Galerkin method and present some numerical experiments to validate our analysis and to demonstrate the effectiveness of our algorithm. The rest of this paper is organized as follows. In the next section, some preliminaries needed in this paper are presented. In §3, the error estimates of approximate eigenvalues and eigenfunctions are analyzed. In §4, we describe the details for an efficient implementation of the algorithm. In §5 we present several numerical experiments to demonstrate the accuracy and efficiency of our method. Finally, in §6 we give some concluding remarks.

2

Preliminaries

Following [7], we consider in this article a particular class of nonlinear eigenvalue problems arising in the study of variational models of the form  v 2 dx = 1}, (2.1) inf{E(v) : v ∈ X, Ω

where X = H01 (Ω) with Ω being a bounded domain, and the energy functional E is of the form  1 1 E(v) = a(v, v) + F (v 2 (x))dx (2.2) 2 2 Ω with

 a(u, v) =

 Ω

(A∇u) · ∇vdx +

Ω

V uvdx.

(2.3)

We make the following assumptions: A(x) is symmetric, and A ∈ (L∞ (Ω))d×d ;

(2.4)

∃ α > 0 s.t. ξ T A(x)ξ ≥ α|ξ|2 for all ξ ∈ Rd and x ∈ Ω;

(2.5)

2

V ∈ L (Ω); 1

(2.6) 2



F ∈ C ([0, +∞), R) ∩ C ((0, +∞), R) and F > 0 on (0, +∞);

(2.7)



(2.8)

∃ 0 ≤ q < 2 and C ∈ R+ s.t. |F (t)| ≤ C(1 + t ) ∀t ≥ 0; 

q



F (t), F (t) bounded in [0, +∞).

(2.9)

In order to simplify the notation, we let f (t) = F  (t), and ω = v 2 . We can then reformulate (2.1) as  √ inf{G(ω) : ω ≥ 0, ω ∈ X, ω = 1}, (2.10) Ω

where 1 √ √ 1 G(ω) = a( ω, ω) + 2 2

2

 Ω

F (ω)dx.

It is shown in [7] that, under assumptions (2.4)-(2.8), (2.10) has a unique solution ω0 and √ ateaux differentiable (2.1) has exactly two solutions: u = ω0 and −u. Moreover, E is Gˆ on X for all v ∈ X, E  (v) = Av v with Av = −div(A∇·) + V + f (v 2 ) being a self-adjoint operator on L2 (Ω). Thus, the function u is a solution of the EulerLagrange equation Au u − λu, v X  ,X = 0, ∀v ∈ X

(2.11)

for some λ ∈ R. We assume that {XN } is a sequence of approximation spaces for X such that lim

min u − vN H 1 = 0.

N →∞ vN ∈XN

Then the discrete variational approximation of (2.1) is:  2 vN = 1}. inf{E(vN ) : vN ∈ XN ,

(2.12)

(2.13)

Ω

Problem (2.13) has at least one minimizer uN , which satisfies AuN uN − λN uN , vN X  ,X = 0, ∀vN ∈ XN ,

(2.14)

for some λN ∈ R.

3

Error estimates

We will establish our main results in this section. We start with a basic error analysis for general approximation spaces, and then derive improved error estimates for the LegendreGalerkin approximation.

3.1

Basic Error Analysis

Let u be the unique positive solution of (2.1) and let uN be a minimizer of the discrete problem (2.13). Since if uN is a minimizer of (2.13), so is −uN , we can assume that (uN , u) ≥ 0. We also introduce the bilinear form E  (u) defined on X × X by   (3.1) E (u)v, w X  ,X = Au v, w X  ,X + 2 f  (u2 )u2 vw. Ω

By using the same arguments as in the proof of Lemma 1 in [7], we can obtain the following Lemma: Lemma 3.1 Under assumptions (2.4)-(2.9) and (2.12), there exist β > 0 and M ∈ R+ such that for all v ∈ X, 0 ≤ (Au − λ)v, v X  ,X ≤ M v2H1 , βv2H1



≤ (E (u) − λ)v, v X  ,X ≤ 3

(3.2) M v2H1 .

(3.3)

And there exists γ > 0 such that for all N > 0, γuN − u2H1 ≤ (Au − λ)(uN − u), (uN − u) X  ,X .

(3.4)

Lemma 3.2 Under assumptions (2.4)-(2.9) and (2.12), it holds that lim uN − uH 1 = 0.

N →∞

(3.5)

Proof. Following [7], we can derive from the definition E(u) that 1 1 E(uN ) − E(u) = Au uN , uN X  ,X − Au u, u X  ,X 2 2 1 + F (u2N ) − F (u2 ) − f (u2 )(u2N − u2 ) 2 Ω 1 = (Au − λ)(uN − u), (uN − u) X  ,X 2 1 + F (u2N ) − F (u2 ) − f (u2 )(u2N − u2 ). 2 Ω From (3.4) and the convexity of F , we have E(uN ) − E(u) ≥

γ uN − u2H 1 . 2

Let πN u ∈ XN be such that u − πN uH 1 = min{u − vN H 1 , vN ∈ XN }. We deduce from (2.12) that πN u converges to u in X when N → ∞. The functional E being strongly continuous on X, then when N → ∞ we have uN − u2H 1 ≤

2 2 (E(uN ) − E(u)) ≤ (E(πN u) − E(u)) → 0. γ γ

Thus, we have lim uN − uH 1 = 0.

N →∞

2

From Lemma 3.2 we know that there exist N1 > 0 such that for all N > N1 , 1 (3.6) uN H 1 ≤ 2uH 1 , uN − uH 1 ≤ . 2 Let Ω = Ω1 ∪ Ω2 , and Ω1 ∩ Ω2 = Ø, such that ∀x ∈ Ω1 , u(uN − u) ≥ 0 and ∀x ∈ Ω2 , u(uN − u) < 0. Since 2 f (u2N ) − f (u2 ) = f  (ξN )(u2N − u2 ), 2 is between u2 and u2 , and f  (t) locally bounded in [0, +∞) and f  (t) > 0 on where ξN N (0, +∞), then from (3.6) there exist non-negative constant αN , βN and M > 0 such that αN < M, βN < M and    2 2 2f (ξN )u (u(uN − u)) = αN u(uN − u), Ω1 Ω1   2 2f  (ξN )u2 (u(uN − u)) = βN u(uN − u). Ω2

Ω2

4

Theorem 1 Under assumptions (2.4)-(2.7), (2.9) and (2.12), it holds |λN − λ| ≤ C(uN − u2H 1 + uN − uL2 ),

(3.7)

uN − uH 1 ≤ C min vN − uH 1 ;

(3.8)

vN ∈XN

In addition, if αN ≤ βN , we have |λN − λ| ≤ CuN − u2H 1 ,

(3.9)

where C is a constant independent of N . Proof. We shall first prove (3.7). Since XN ⊂ H01 (Ω), we derive from (2.1), (2.11), (2.13) and (2.14) that λN ≥ λ. On the other hand, by direct calculation, we have λN − λ = AuN uN , uN X  ,X − Au u, u X  ,X   =a(uN , uN ) − a(u, u) + f (u2N )u2N − f (u2 )u2 Ω Ω  =a(uN − u, uN − u) + 2a(u, uN − u) + f (u2N )u2N − f (u2 )u2 Ω Ω   =a(uN − u, uN − u) + 2λ u(uN − u) − 2 f (u2 )u(uN − u) Ω Ω   + f (u2N )u2N − f (u2 )u2 Ω

Ω

=a(uN − u, uN − u) − λuN − u2L2    − 2 f (u2 )u(uN − u) + f (u2N )u2N − f (u2 )u2 Ω Ω Ω  = (Au − λ)(uN − u), (uN − u) X  ,X + u2N (f (u2N ) − f (u2 )). Ω

We estimate below the two terms in the last line. Since u2N (u2N − u2 ) = (2u2 + 2uuN + u2N )(uN − u)2 + 2u3 (uN − u), we have

 Ω

u2N (f (u2N ) − f (u2 ))dx =

 Ω 

2 f  (ξN )(2u2 + 2uuN + u2N )(uN − u)2 dx

+2 Ω

2 f  (ξN )u3 (uN − u)dx.

The above two terms can be estimated as follows:  2 f  (ξN )(2u2 + 2uuN + u2N )(uN − u)2 dx Ω

≤ C(u2 + u2N )(uN − u)L2 uN − uL2 ≤ C(u2 (uN − u)L2 + u2N (uN − u)L2 )uN − uL2 ≤ C(u2L6 uN − uL6 + uN 2L6 uN − uL6 )uN − uL2

≤ C(c36 u2H 1 uN − uH 1 + 4c36 u2H 1 uN − uH 1 )uN − uL2 ≤ CuN − uH 1 uN − uL2 ≤ CuN − u2H 1 , 5

 2 Ω

2 f  (ξN )u3 (uN − u)dx ≤ Cu3L6 uN − uL2

≤ C(c36 u3H 1 )uN − uL2 ≤ CuN − uL2 . Therefore, we have  Ω

u2N (f (u2N ) − f (u2 ))dx ≤ C(uN − u2H 1 + uN − uL2 ),

where c6 is the Soblev constant in vL6 ≤ c6 vH 1 , vN

Next, we will evaluate the ∈ XN ,

H 1 -norm

∀v ∈ X. Therefore, we obtain (3.7).

of the error uN − u. We first notice that for all

uN − uH 1 ≤ uN − vN H 1 + vN − uH 1 .

(3.10)

From (3.3) of Lemma 3.1 we have uN − vN 2H 1 ≤ β −1 (E  (u) − λ)(uN − vN ), (uN − vN ) X  ,X = β −1 (E  (u) − λ)(uN − u), (uN − vN ) X  ,X

+ β −1 (E  (u) − λ)(u − vN ), (uN − vN ) X  ,X . We proceed below in three steps. Step 1: Estimation of (E  (u) − λ)(uN − u), (uN − vN ) X  ,X . Since (E  (u) − λ)(uN − u), (uN − vN ) X  ,X  = − (f (u2N )uN − f (u2 )uN − 2f  (u2 )u2 (uN − u))(uN − vN )dx Ω  + (λN − λ) uN (uN − vN )dx, 

Ω

2 2  2 2 we  only need to estimate Ω (f (uN )uN − f (u )uN − 2f (u )u (uN − u))(uN − vN )dx and Ω uN (uN − vN )dx. For all vN ∈ XN such that vN L2 = 1, we have   1 uN (uN − vN )dx = 1 − uN vN dx = uN − vN 2L2 . 2 Ω Ω

6

In addition, we have  | (f (u2N )uN − f (u2 )uN − 2f  (u2 )u2 (uN − u))(uN − vN )dx| Ω 2 ≤ |f  (ξN )(uN + u)uN − 2f  (u2 )u2 )(uN − u)| · |uN − vN |dx Ω  2 2 2 ≤ |(f  (ξN )(uN + u)uN − 2f  (ξN )u2 + 2f  (ξN )u2 Ω

 ≤ ≤ ≤ ≤

− 2f  (u2 )u2 )(uN − u)| · |uN − vN |dx

Ω Ω Ω Ω

2 2 |(f  (ξN )(uN − u)(uN + 2u) − 2u2 (f  (ξN ) − f  (u2 )))(uN − u)| · |uN − vN |dx 2 2 2 |(f  (ξN )(uN − u)(uN + 2u) − 2u2 f  (ξ1N )(ξN − u2 ))(uN − u)| · |uN − vN |dx 2 2 2 (|f  (ξN )(uN + 2u)(uN − u)| + |2u2 f  (ξ1N )(ξN − u2 )|) · |uN − u| · |uN − vN |dx 2 2 (|f  (ξN )(uN + 2u)(uN − u)| + |2u2 f  (ξ1N )(u2N − u2 )|) · |uN − u| · |uN − vN |dx

2 2 (|f  (ξN )(uN + 2u)| + |2u2 f  (ξ1N )(uN + u)|) · (uN − u)2 · |uN − vN |dx Ω ≤ C (|uN + 2u| + |u2 (uN + u)|) · (uN − u)2 · |uN − vN |dx,

≤

Ω

2 is between ξ 2 and u2 . Since for all N > N and all v ∈ X , where ξ1N 1 N N N  |uN + 2u| · (uN − u)2 · |uN − vN |dx Ω

≤ uN + 2uL2 (uN − u)2 (uN − vN )L2 ≤ (uN L2 + 2uL2 )uN − u2L6 uN − vN L6

≤ (CuN L6 + 2uL2 )uN − u2L6 uN − vN L6

≤ (Cc6 uN H 1 + 2uL2 )c36 uN − u2H 1 uN − vN H 1 ≤ (2Cc6 uH 1 + 2uL2 )c36 uN − u2H 1 uN − vN H 1 ≤ CuN − u2H 1 uN − vN H 1 , and

 Ω

|u2 (uN + u)| · (uN − u)2 · |uN − vN |dx ≤ u2 (uN + u)L2 (uN − u)2 (uN − vN )L2 ≤ (u2 uN L2 + u3 L2 )(uN − u)2 (uN − vN )L2 ≤ (u2L6 uN L6 + u3L6 )uN − u2L6 uN − vN L6

≤ c36 (u2H 1 uN H 1 + u3H 1 )uN − u2L6 uN − vN L6 ≤ c66 (2u3H 1 + u3H 1 )uN − u2H 1 uN − vN H 1 ≤ CuN − u2H 1 uN − vN H 1 ,

7

we derive from (3.7) that | (E  (u) − λ)(uN − u), (uN − vN ) X  ,X | ≤ C(uN − u2H 1 uN − vN H 1 + (uN − u2H 1 + uN − uL2 )uN − vN 2H 1 ). Step 2: Estimation of (E  (u) − λ)(u − vN ), (uN − vN ) X  ,X . By direct calculations, we find | Au (u − vN ), uN − vN X  ,X |  ≤ | (−div(A∇(u − vN )) + V (u − vN ) + f (u2 )(u − vN ))(uN − vN )dx| Ω ≤ | (A∇(u − vN )∇(uN − vN ) + V (u − vN )(uN − vN ) + f (u2 )(u − vN ))(uN − vN )dx| Ω  ≤ AL∞ ∇(u − vN )L2 ∇(uN − vN )L2 + |V (u − vN )(uN − vN )|dx Ω  + |f (u2 )(u − vN )(uN − vN )|dx Ω

≤ AL∞ u − vN H 1 uN − vN H 1 + V L2 (u − vN )(uN − vN )L2 + f (u2 )L∞ u − vN L2 uN − vN L2 ≤ AL∞ u − vN H 1 uN − vN H 1 + CV L2 (u − vN )(uN − vN )L3 + f (u2 )L∞ u − vN L2 uN − vN L2 ≤ AL∞ u − vN H 1 uN − vN H 1 + CV L2 u − vN L6 uN − vN L6 + f (u2 )L∞ u − vN L2 uN − vN L2 ≤ AL∞ u − vN H 1 uN − vN H 1 + Cc26 V L2 u − vN H 1 uN − vN H 1 + f (u2 )L∞ u − vN L2 uN − vN L2 ≤ (AL∞ + Cc26 V L2 + f (u2 )L∞ )u − vN H 1 uN − vN H 1 ≤ Cu − vN H 1 uN − vN H 1 , and

 2

 |f  (u2 )u2 (u − vN )(uN − vN )|dx + λ |(u − vN )(uN − vN )|dx Ω Ω  ≤C |u2 (u − vN )(uN − vN )|dx + λu − vN L2 uN − vN L2 Ω  3 3 2 2 ≤ CuL6 ( (u − vN ) 2 (uN − vN ) 2 dx) 3 + λu − vN L2 uN − vN L2 ≤ ≤ ≤

Ω 2 2 Cc6 uH 1 u − vN L2 uN − vN L6 + λu − vN L2 uN − vN L2 Cc36 u2H 1 u − vN L2 uN − vN H 1 + λu − vN H 1 uN − vN H 1 (Cc36 u2H 1 + λ)u − vN H 1 uN − vN H 1

≤ Cu − vN H 1 uN − vN H 1 .

8

Therefore, we have | (E  (u) − λ)(u − vN ), (uN − vN ) X  ,X |  ≤ | (Au (u − vN ), uN − vN X  ,X | + 2 |f  (u2 )u2 (u − vN )(uN − vN )|dx Ω  + λ |(u − vN )(uN − vN )|dx ≤ Cu − vN H 1 uN − vN H 1 . Ω

Step 3: Estimation of uN − uH 1 . From Step 1 and Step 2, we have uN − vN H 1 ≤ C(uN − u2H 1 + u − vN H 1

+ (uN − u2H 1 + uN − uL2 )uN − vN H 1 ).

We derive from (3.5) and (3.6) that there exist N2 > N1 such that for all N ≥ N2 , C(uN − u2H 1 + uN − uL2 ) ≤ γ <

1 2

1 CuN − uH 1 ≤ γ < , 2

(3.11) (3.12)

where γ is a constant independent of N . Thus, from (3.10) for all N > N2 we have uN − uH 1 ≤ C

min

vN ∈XN ,vN L2 =1

vN − uH 1 .

(3.13)

We let u0N be a minimizer of the following minimization problem min vN − uH 1 .

vN ∈XN

We know from (2.12) that u0N converges to u in H 1 when N → ∞. In addition, we have min

vN ∈XN ,vN L2 =1

vN − uH 1 ≤ 

u0N − uH 1 u0N L2

≤ u0N − uH 1 + ≤ u0N − uH 1 + ≤ (1 + ≤ (1 +

u0N H 1 |1 − u0N L2 | u0N L2 u0N H 1 0 u − uL2 u0N L2 N

u0N H 1 )u0N − uH 1 u0N L2

u0N H 1 ) min vN − uH 1 . u0N L2 vN ∈XN

For N > N2 > N1 , we have 1 u0N − uH 1 ≤ uN − uH 1 ≤ , 2 1 u0N − uL2 ≤ uN − uH 1 ≤ . 2 9

Then we have u0N H 1 ≤ u0N − uH 1 + uH 1 ≤

1 + uH 1 . 2

Since 1 = uL2 ≤ u0N − uL2 + u0N L2 ≤

1 + u0N L2 , 2

then we have 1 u0N L2 ≥ . 2 Then we can get 1+

u0N H 1 ≤ 2(uH 1 + 1). u0N L2

Thus, (3.8) is proved. Finally, if αN ≤ βN , then we have     2 2  2 2 2 2f (ξN )u (u(uN − u)) = 2f (ξN )u (u(uN − u)) + 2f  (ξN )u2 (u(uN − u)) Ω Ω1 Ω2   u(uN − u) + βN u(uN − u) = αN Ω1 Ω2   u(uN − u) + u(uN − u)) ≤ βN ( Ω1 Ω2  1 u(uN − u) = − βN uN − u2L2 . = βN 2 Ω Then from (3.2) in Lemma 3.1, we obtain (3.9). 2 Remark 3.1. The optimal eigenvalue error estimate (3.9) is proved without additional smoothness assumption but with the assumption αN ≤ βN which is not easy to verify. However, the numerical results in §5 indicates that this optimal eigenvalue error estimate holds, at least for the tested cases.

3.2

Error estimates for Legendre-Galerkin method

The error estimates in Theorem 1 is proved under very general assumptions on the approximation space. In this subsection, we shall fix Ω = I d (d = 1, 2, 3) with I = (−1, 1), and consider the Legendre-Galerkin approximation with XN = PN (I d ) ∩ H01 (I d ) where PN stand for the set of all polynomials of at most degree N in each direction. Then (2.12) is obviously satisfied. 1,0 : H01 (I d ) → XN by We define the projection operator πN 1,0 (∇(πN u − u), ∇v) = 0,

10

v ∈ XN ,

and recall the following result (cf. Remark 2.16 of [5]): Lemma 3.3 Let r ∈ N0 with r ≥ 1. If f ∈ H01 (I d ) ∩ H r (I d ), then, for N ≥ 1, 1,0 f H 1 (I d ) ≤ CN −r+1 f H r (I d ) , f − πN

where C is a constant independent of N . Combining the above results with those in Theorem 1, we derive immediately the following: Theorem 2 Under assumptions (2.4)-(2.7) and (2.9), if u ∈ H01 (I d ) ∩ H m (I d ), then for m ≥ 1, N ≥ 1, we have |λN − λ| ≤ CN −m+1 uH m , uN − uH 1 ≤ CN

−m+1

(3.14)

uH m ;

(3.15)

In addition, if αN ≤ βN , we have |λN − λ| ≤ CN 2(−m+1) (u2H m + uH m ),

(3.16)

where C is a constant independent of N .

4

Implementation details and numerical results

In this section, we present an efficient implementation of the Legendre-Galerkin method for (2.14), and present some numerical experiments to validate our error analysis.

4.1

Implementation of the Legendre-Galerkin method

From equation (2.11) and the constraint u2L2 = 1, we obtain an equivalent nonlinear eigenvalue problem − div(A∇u) + V u + f (u2 )u = λu, u2L2 = 1 in Ω,

(4.1)

u|∂Ω = 0.

(4.2)

The weak form of (4.1) and (4.2) is: Find (u, λ) ∈ X × R such that (A∇u, ∇v) + (V u, v) + (f (u2 )u, v) = λ(u, v), u2L2

∀v ∈ X,

= 1.

(4.3) (4.4)

Then the discrete form of (4.3) and (4.4) is: Find (uN , λN ) ∈ XN × R such that (A∇uN , ∇vN ) + (V uN , vN ) + (f (u2N )uN , vN ) = λN (uN , vN ), uN 2L2

∀vN ∈ XN ,

= 1.

(4.5) (4.6)

To solve the above nonlinear eigenvalue problem, we use the Picard iteration as follows: p p 2 p (A∇upN , ∇vN ) + (V upN , vN ) + (f ((up−1 N ) )uN , vN ) = λN (uN , vN ), ∀vN ∈ XN ,

upN 2L2

= 1,

(4.7) (4.8)

11

with initial guess determined by (A∇u0N , ∇vN ) + (V u0N , vN ) = λ0N (u0N , vN ), ∀vN ∈ XN , u0N 2L2

= 1.

(4.9) (4.10)

To simplify the presentation, we shall only consider the two-dimensional case although higher-dimensional case can be dealt with similarly. Let ϕk = Lk (x) − Lk+2 (x) (k = 0, 1, · · · , N − 2), where Lk (x) denotes the Legendre polynomial of degree k. Then, we have XN = span{ϕi (x)ϕj (y), i, j = 0, 1, · · · , N − 2} [15]. We write upN

=

N −2 

upij ϕi (x)ϕj (y).

(4.11)

i,j=0

Then, we can reduce (4.7)-(4.8) and (4.9)-(4.10) to generalized eigenvalue problems: up = λpN B¯ up , p = 1, 2, · · · , (S + M(V ) + M(¯ up−1 ))¯

(4.12)

u0 , (S + M(V ))¯ u0 = λ0N B¯

(4.13)

and

respectively, where the stiff matrix S and mass matrix B are sparse [15]. The matrices M(V ) v and M(¯ up−1 )¯ v and M(¯ up−1 ) are in general full but their matrix-vector product (M(V ))¯ can be efficiently computed by using a pseudo-spectral approach. Since only a few smallest eigenvalues are mostly interesting in real applications, it is most efficient to solve (4.12) (resp. (4.13)) using iterative eigen solvers such as shifted inverse power method (cf., for instance, [19]) which requires solving, repeatedly for different righthand side f (resp. f˜), ((S + M(V ) + M(¯ up−1 )) − λap B)up = f (resp. ((S + M(V )) − λa0 B)u0 = f˜),

(4.14)

where λap (resp. λa0 ) is some approximate value for the eigenvalue λpN (resp. λ0N ). The above system can be efficiently solved by the Schur-complement approach, we refer to [20] for a detailed description on a related problem. In summary, the approximate eigenvalue problem (4.12) (resp. (4.13)) can be solved very efficiently.

4.2

Numerical experiments

We now perform some numerical tests to compute eigenvalues and eigenfunctions of (3.14)(3.16). All numerical tests are performed using MATLAB 2015b. 4.2.1

One dimensional case

¯ = [−1, 1] as our first example. Numerical We take A = 12 , V (x) = 12 x2 , f (u2 ) = |u|2 and Ω results for the first eigenvalue with different N and iteration step L are listed in Tables 5.1. Since the exact eigen-pairs are unknown, we computed the reference solution, λL N,1 and the L associated numerical eigenfunction uN,1 , with N = 40 and iteration step L = 30.

12

Table 5.1. Numerical approximation to λ1 for different N and L in 1D.

N 10 15 20 25

L = 10 2.038315754257193 2.038315716117307 2.038315716123191 2.038315716123191

L = 15 2.038315754334226 2.038315716190542 2.038315716196435 2.038315716196430

L = 20 2.038315754334225 2.038315716190543 2.038315716196433 2.038315716196430

L = 25 2.038315754334222 2.038315716190548 2.038315716196430 2.038315716196434

We plot the error graphs of numerical approximations to the first eigenvalues λ1 and the associated eigenfunction u1 with different N and L in Figure 1 and Figure 3, respectively. 2 L In addition, in order to compare the convergence rates of u1 − uL N,1 H 1 and u1 − uN,1 L2 , 30 L 30 L 2 30 we also plot the error graphs of u40,1 − uN,1 L2 and u40,1 − uN,1 H 1 /u40,1 − uL N,1 L2 for the first eigenvalue in Figure 2, 4. We see from Table 5.1 and Figure 1 that numerical eigenvalues achieve fifteen-digit accuracy with N ≥ 20 and L ≥ 15. We know from Figure 1 that the numerical results are in agreement with the theoretical results, i.e., achieve 2 spectral accuracy. From Figure 4, the convergence rate of u1 − uL N,1 H 1 is higher than that of u1 − uL N,1 L2 . However, we know from Figure 1 and Figure 3 that the convergence 2 − λ1 | and uL rates of |λL N,1 N,1 − u1 H 1 are almost the same, which show that correctness of optimal error estimation in Theorem 1. −2

10-2

10 L=15 L=20 L=25

10-4

L=15 L=20 L=25

−4

10

−6

10-6 N,1 L

||u30 −uL ||

10-8

40,1

Error

2

10

10-10

−8

10

−10

10

−12

10-12

10

10-14

10

10-16

−14

−16

10

6

8

10

12

14

16

18

20

5

10

15

20

25

30

N

N

L Figure 1: Errors between numerical solutions Figure 2: The error figure of u30 40,1 − uN,1 L2 for different N and L in 1D. and the reference solution for λ1 in 1D.

−2

0

10

10 L=15 L=20 L=25

−4

10

−6

||u30 −uL ||2 1 /||u30 −uL ||

N,1 L

2

10

40,1

−10

N,1 H

10

40,1

N,1 H

||u30 −uL ||2 1

−8

10

−12

40,1

10

−14

10

−4

10

−6

10

−8

10

−10

10

−12

10

−16

10

−18

10

L=15 L=20 L=25

−2

10

−14

6

8

10

12

14

16

18

10

20

N

5

10

15

20

25

30

N

L 2 Figure 3: The error figure of u30 The error figure of u30 40,1 − uN,1 H 1 Figure 4: 40,1 − L 2 30 for different N and Lin 1D. uN,1 H 1 /u40,1 − uL  for different N and 2 N,1 L L in 1D.

Remark 4.1. For the non-periodic case, the numerical method in reference [7] is mainly 13

based on the finite element discretization. The numerical method in this paper is based on spectral Galerkin approximation, thus, when the solution is smooth enough, the numerical solutions have spectral accuracy. In addition, we can observe from the numerical results in this paper that the numerical solutions have spectral accuracy. Compared with the numerical results in reference [7] , the accuracy of the numerical solution in this paper is much higher in the same degree of freedom. 4.2.2

Two dimensional case

¯ = [−1, 1]2 as our second We take A = 12 E, V (x1 , x2 ) = 12 (x21 + x22 ) , f (u2 ) = |u|2 and Ω example, where E is the identity matrix. Numerical results for the first eigenvalue with different N and iteration step L are listed in Tables 5.2. Since the exact eigen-pairs are unknown, we computed the reference solution, λL N,1 and the associated numerical eigenL function uN,1 , with N = 40 and iteration step L = 30. Table 5.2. Numerical approximation to λ1 for different N and L in 2D.

N 10 15 20 25 30

L=6 3.150223951725874 3.150223942247921 3.150223942247918 3.150223942247905 3.150223942247945

L = 12 3.150224069249848 3.150224059885169 3.150224059885197 3.150224059885206 3.150224059885078

L = 18 3.150224069249997 3.150224059885314 3.150224059885376 3.150224059885377 3.150224059885060

L = 24 3.150224069249999 3.150224059885318 3.150224059885364 3.150224059885364 3.150224059885223

We also plot the error graphs of numerical approximations to the first eigenvalue λ1 and the associated eigenfunction u1 with different N and L in Figure 5 and Figure 6, respectively. We see from Table 5.2 and Figure 5 that numerical eigenvalues achieve at least thirteen-digit accuracy with N ≥ 15 and L ≥ 12. We know from Figure 5 that the numerical results are in agreement with the theoretical results, i.e., achieve spectral accuracy. However, we know L 2 from Figure 5 and Figure 6 that the convergence rates of |λL N,1 − λ1 | and uN,1 − u1 H 1 are almost the same, which also show that correctness of optimal error estimation in Theorem 1. In order to further demonstrate the convergence of the approximate eigenfunctions, we plot the graphs of the first eigenfunction with N = 15 and L = 18 in Figure 7 and the graphs of the first reference eigenfunction with N = 40 and L = 30 in Figure 8, respectively.

−4

0

10

10 L=12 L=18 L=24

L=12 L=18 L=24

−6

10

−5

10

||u30 −uL ||2 1

N,1 H

−10

10

40,1

40,1

N,1

|λ30 −λL |

−8

10

−10

10

−15

10

−12

10

−14

10

−20

6

8

10

12

14

16

18

10

20

N

6

8

10

12

14

16

18

20

N

L 2 Figure 5: Errors between numerical solutions Figure 6: The error figure of u30 40,1 − uN,1 H 1 for different N and L in 2D. and the reference solution for λ1 in 2D.

14

0.6

0.6 u30

u18

40,1

1 0.8

15,1

1 0.8

0.4 0.2

0.4 0.2

0 1

0 1 0.5

0.5

1 0.5

0

y

1

−1

0

−0.5

−0.5 −1

0.5

0

0

−0.5

y

x

−0.5 −1

−1

x

Figure 7: The graphs of the first eigenfunction Figure 8: The graphs of the first eigenfunction with N = 15 and L = 18 in 2D. with N = 40 and L = 30 in 2D.

4.3

Summary

We considered numerical approximations and error estimates for a nonlinear elliptic eigenvalue problem. Spectral accuracy error bounds are established for numerical eigenvalues and eigenfunctions. Numerical tests demonstrate that the method achieves high accuracy with relatively small number of unknowns. To simplify the analysis, we have restricted our analysis to rectangle and cubic domains. However, the approach presented in this paper can be extended to more general domains by using spectral-element method. Acknowledgements. This work was supported in part by the National Natural Science Foundation of China (No. 11661022, 11471031, 91430216), NASF U1530401, and the US National Science Foundation grant DMS-1419040.

References [1] Babuˇ ska I., Osborn J., Eigenvalue problems. Handbook of numerical analysis, 1991, 2: 641-787. [2] Bao W. and Du Q., Computing the ground state solution of Bose-Einstein condensates by a normalized gradient flow. SIAM Journal on Scientific Computing, 2004, 25(5): 1674-1697. [3] Bao W., Jaksch D. and Markowich P. A., Numerical solution of the Gross-Pitaevskii equation for Bose Einstein condensation. Journal of Computational Physics, 2003, 187(1): 318-342. [4] Bao W., Chern I. L. and Lim F. Y., Efficient and spectrally accurate numerical methods for computing ground and first excited states in Bose-Einstein condensates. Journal of Computational Physics, 2006, 219(2): 836-854. [5] Bernardi C. and Maday Y., Approximations spectrales de problemes aux limites elliptiques. Springer, 1992. [6] Cances E., Defranceschi M., Kutzelnigg W. et al., Computational quantum chemistry: a primer. Handbook of numerical analysis, 2003, 10: 3-270.

15

[7] Cances E., Chakir R. and Maday Y. Numerical analysis of nonlinear eigenvalue problems. Journal of Scientific Computing, 2010, 45(1-3): 90-117. [8] Cances E., Chakir R., and Maday Y. Numerical analysis of the planewave discretization of some orbital-free and Kohn-Sham models, M2AN, 2012 46: 341-388. [9] Chen H., Gong X. and Zhou A., Numerical approximations of a nonlinear eigenvalue problem and applications to a density functional model. Mathematical Methods in the Applied Sciences, 2010, 33(14): 1723-1742. [10] Chen H., He L. and Zhou A., Finite element approximations of nonlinear eigenvalue problems in quantum physics. Computer Methods in Applied Mechanics and Engineering, 2011, 200(21): 1846-1865. [11] Chen H., Gong X., He L. et al., Numerical analysis of finite dimensional approximations of Kohn-Sham models. Advances in Computational Mathematics, 2013: 1-32. [12] Neese F., Prediction of electron paramagnetic resonance g values using coupled perturbed Hartree-Fock and Kohn-Sham theory. The Journal of Chemical Physics, 2001, 115(24): 11080-11096. [13] Motamarr P., Nowak M. R., Leiter K., Knap J., Gavini V., Higher-order adaptive finiteelement methods for Kohn-Sham density functional theory. Journal of Computational Physics, 2013, 253(C):308-343 [14] Pitaevskii L. P. and Stringari S.; Bose-einstein condensation. Oxford University Press, 2003. [15] Shen J., Efficient spectral-Galerkin method I. Direct solvers of second-and fourth-order equations using Legendre polynomials. SIAM Journal on Scientific Computing, 1994, 15(6): 1489-1505. [16] Tkatchenko A. and Scheffler M., Accurate molecular van der Waals interactions from ground-state electron density and free-atom reference data. Physical Review Letters, 2009, 102(7): 073005. [17] Zhou A., An analysis of finite dimensional approximations for the ground state solution of Bose-Einstein condensates. Nonlinearity, 2003, 17(2): 541. [18] Zhou A., Finite dimensional approximations for the electronic ground state solution of a molecular system. Mathematical methods in the applied sciences, 2007, 30(4): 429-447. [19] Golub G. H. and Van Load C. F., Matrix Computations. The John Hopkins University Press, Baltimore, 1989. [20] Kwan Y. and Shen J., An efficient direct parallel elliptic solver by the spectral element method. J. Comput. Phys., 2007, 225: 1721-1735.

16