Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Convergence of a normalized gradient algorithm for computing ground states

Convergence of a normalized gradient algorithm for computing ground states Abstract We consider the approximation of the ground state of the one-dimensional cubic nonlinear Schrödinger equation by a normalized gradient algorithm combined with a linearly implicit time integrator, and a finite difference space approximation. We show that this method, also called the imaginary time evolution method in the physics literature, is locally convergent, and we provide error estimates: for initial data in a neighborhood of the ground state, the algorithm converges exponentially toward a modified soliton that is a space discretization of the exact soliton, with error estimates depending on the discretization parameters. 1. Introduction The goal of this article is to give a convergence proof for a normalized gradient algorithm used to numerically compute ground states of Schrödinger equations fulfilling symmetry and coercivity conditions as considered in the seminal works of Weinstein (1985) and Grillakis et al. (1987, 1990). This algorithm is also called the imaginary time method in the physics literature; see for instance Edwards & Burnett (1995), Chiofalo et al. (2000), Adhikari (2000a,b), Bao & Tang (2003) and the references therein. Let us describe the algorithm in the case of the focusing cubic nonlinear Schrödinger equation   i∂tψ=−12Δψ−|ψ|2ψ, (NLS) set on $${\mathbb{R}}$$, where $$\psi(t,x)$$ depends on space variables $$x \in {\mathbb{R}}$$. With this equation is associated the energy   H(ψ,ψ¯)=14∫R|∇ψ|2−|ψ|4, (1.1) which is preserved by the flow of (NLS) for all times. The equation (NLS) can be written as   i∂tψ=−12Δψ−|ψ|2ψ=2∂H∂ψ¯(ψ,ψ¯). In the rest of this article, the notation $$\nabla H$$ will denote the $$L^2$$ derivative of the energy $$H$$ with respect to real functions $$\psi$$. Note that we have for a real function $$u$$,   ∇H(u)=2∂H∂ψ¯(u,u)=−12Δu−u3, the left-hand side denoting the Fréchet derivative of $$H(u)$$ considered as a functional acting on real functions. Note that naturally, $$\nabla H \in H^{-1}$$ with the embedding $$H^{-1} \subset L^2 \subset H^1$$. With this notation, the ground state $$\eta(x)$$ is defined as the unique real symmetric minimizer in $$H^1$$ of the problem   min‖ψ‖L2=1H(ψ). (1.2) In the one-dimensional cubic case considered in this article, explicit computations show that   η(x):=12sech(x2). We denote by $$\lambda$$ the Lagrange multiplier associated with this minimization problem, such that   ∇H(η)=−12Δη−η3=−λη. (1.3) With the ground state $$\eta$$ is associated the solution $$\eta(t,x) = \eta(x) e^{{\rm i} \lambda t}$$ of the time-dependent equation (NLS). Using translation and scaling, the ground state gives rise to a family of explicit solutions of (NLS) that have the property of being orbitally stable in $$H^1$$ (see Weinstein, 1985; Grillakis et al., 1987, 1990; Bambusi et al., 2013). For instance, any solution starting close to $$\eta(x)$$ will remain close to the manifold $$\{\,{\it e}^{{\rm i}\alpha} \eta ( x - c)\, | \alpha,c \in {\mathbb{R}}\} \subset H^1$$, the rotation and translation being natural invariant group actions of the nonlinear Schrödinger equation (NLS). In more general situations, $$\eta$$ is not explicitly known, and one has to rely on numerical simulations to compute it. To this aim, the imaginary time method, which is a nonlinear version of the normalized gradient algorithm, is widely used. The goal of the present article is to analyse the efficiency of this method in the simple case described above (the results obtained are in fact valid in more general situations—essentially all situations where the Grillakis–Shatah–Strauss arguments apply). The time-discretized algorithm consists in defining a sequence a functions $$\{\psi_n\}_{n \in {\mathbb{N}}}$$ as follows: (i) For a given $$\psi_n$$, an intermediate function $$\psi_{n}^*$$ is defined as a numerical approximation of the solution of the parabolic equation   ∂tψ=12Δψ+|ψ|2ψ=−∇H(ψ), (1.4) starting in $$\psi_n$$, over a time interval $$[0,\tau]$$, where $$\tau$$ is a given time step. To compute $$\psi_n^*$$, we will see that the best results are given by the linearly implicit method   ψn∗=ψn−τ∇H^(ψn,ψn∗), (1.5) where   −∇H^(ψn,ψn∗)=12Δψn∗+|ψn|2ψn∗. (ii) Then, we define the normalized function   ψn+1=ψn∗‖ψn∗‖L2. (1.6) Note that the algorithm presented above preserves the real nature of the function $$\psi$$: if $$\psi_0$$ is real valued then for all $$n \in {\mathbb{N}}$$, $$\psi_n$$ is real, and the same holds true for symmetric functions satisfying $$\psi_n(-x) = \psi_n(x)$$. Our results can be summarized as follows: In the semidiscrete case described above, if the initial data $$\psi_0$$ is real, symmetric and sufficiently close to $$\eta$$ then $$\psi_n$$ converges exponentially toward $$\eta$$ in $$H^1$$. In the fully discrete case, where the discretization in space is made by finite differences in space with mesh $$h$$ combined with a Dirichlet cut-off for large values of $$x$$ (namely $$x \leq Kh$$), then we prove exponential convergence toward a modified soliton $$\eta_{h,K}$$ that is $$\mathcal{O}( h + \frac{1}{h^2} {\it \,{e}}^{- C_1 Kh})$$ close to the exact soliton $$\eta$$. The main property explaining the excellent performance of this method is the fact that the linearly implicit numerical scheme exactly preserves the ground state: if $$\psi_n = \eta$$ then $$\psi_{n+1} = \eta$$, and the same holds true for discrete-in-space ground states satisfying a discrete version of (1.3). This fact is very general: it holds for any linearly implicit scheme applied to a semilinear partial differential equation (PDE). This important feature of course makes this scheme much more attractive than other possible schemes, where the nonlinearity would be approached by   −∇H^(ψn,ψn∗)=12Δψn∗+{|ψn|2ψn(semi−explicit),|ψn∗|2ψn∗(fullyimplicit),  (1.7) despite the fact that fully implicit schemes enjoy the energy-diminishing property for a standard gradient system (see Hairer & Lubich, 2014). As we will discuss later, these schemes still converge, but toward modified ground states $$\eta_\tau = \eta + \mathcal{O}(\tau)$$. As these results would be weaker, with longer detailed proofs, we give only the main arguments to obtain them, keeping a full detailed convergence proof for the linearly implicit scheme. 2. The Hamiltonian near the ground state We work in the space $$V$$ of real symmetric functions of $$H^1$$:   V:={φ∈H1(R,R)|φ(−x)=φ(x)}. We consider the usual $$L^2$$ and $$H^1$$ norms and denote by $$\left\langle{\cdot},{\cdot}\right\rangle$$ the canonical real scalar product on $$L^2$$:   ⟨φ,ψ⟩=∫Rφ(x)ψ(x)ds,‖φ‖L22=⟨φ,φ⟩,‖φ‖H12=‖φ‖L22+‖∂xφ‖L22. 2.1 Coordinates in the neighborhood of the ground state We introduce the set of the functions $$R$$-close to $$\eta$$,   U(R):={φ∈V|‖φ−η‖H1<R}, the set   W:={u∈V|⟨u,η⟩=0} and the map $$\chi$$  χ:R×W→V,(r,u)↦(1+r)η+u. The map $$\chi$$ allows us to use $$(r,u)$$ as coordinates in $$\,\mathcal{U}(R)$$. Observe that $$\chi$$ is smooth with bounded derivatives and so is the inverse $$\chi^{-1}$$, from the explicit formula   χ−1:V→R×W,ψ↦(r(ψ),u(ψ))=(⟨ψ,η⟩−1,ψ−⟨ψ,η⟩η). We will also use the following notation: we define the $$L^2$$ projectors   Pηu=⟨u,η⟩ηandPWu=I−Pη, and we define the function   H~(r,u)=(H∘χ)(r,u). (2.1) We can verify the following relations:   ∂rH~(r,u)=⟨∇H(χ(r,u)),η⟩=η−1(Pη∇H)(χ(r,u))∈R (2.2) and   ∇uH~(r,u)=PW∇H(χ(r,u))∈H−1. (2.3) Before collecting some expressions of the Hamiltonian and the gradient flow in coordinates $$(r,u)$$, let us introduce the following notation. Definition 2.1 Let $$p \geq 1$$. We say that $$R(u) = \mathcal{O}(\left\|{u}\right\|_{H^1}^p)$$ if $$R: H^1 \to H^1$$ is a smooth function such that for all $$B$$, there exists a constant $$C$$ such that for all $$u\in H^1$$, $$\left\|{u}\right\|_{H^1} \leq B$$, we have   ‖R(u)‖H1≤C‖u‖H1p. We will also use the notation $$R(u,v) = \mathcal{O}(\left\|{u,v}\right\|_{H^1}^p)$$ for a smooth function $$R: H^1 \times H^1 \to H^1$$ if for all $$B$$, there exists $$C$$ such that for all $$u,v \in H^1$$ satisfying $$\left\|{u}\right\|_{H^1} \leq B$$ and $$\left\|{v}\right\|_{H^1} \leq B$$ then we have   ‖R(u,v)‖H1≤C(‖u‖H1p+‖v‖H1p). Finally, a function $$u = \mathcal{O}(r)$$ if $$\left\|{u}\right\|_{H^1} \leq C r$$ for $$r$$ small enough and a constant $$C$$ independent of $$u$$. With this notation and the fact that $$\psi = (1 + r) \eta + u$$, we compute using (1.3) that   −∇H(ψ) =12Δψ+ψ3 = (1+r)(−η3+λη)+12Δu+u3 +(1+r)3η3+3(1+r)2η2u+3(1+r)ηu2 =λη+12Δu+3η2u+O(r+‖u‖H12). Note that to obtain the bound, we have used the fact that $$\eta \in H^1$$, as well as the estimate   ‖uv‖H1≤C‖u‖H1‖v‖H1 (2.4) for two functions $$u$$ and $$v$$. We deduce using (2.2) and (2.3) that   ∂rH~(r,u) = −λ−12⟨u,Δη⟩−3⟨η3,u⟩+O(r+‖u‖H12) = −λ−2⟨η3,u⟩+O(r+‖u‖H12), as $$\left\langle{u},{\eta}\right\rangle = 0$$ and $$\left\langle{\eta},{\eta}\right\rangle = 1$$ and   ∇uH~(r,u) = −PW(λη+12Δu+3η2u)+O(r+‖u‖H12) = −12Δu−3η2u+2⟨η3,u⟩η+O(r+‖u‖H12). 2.2 Projection onto the unit $$L^2$$ sphere Let us now define the function $$u\mapsto r(u)$$ from $$W$$ to $${\mathbb{R}}$$ by the implicit relation   ‖χ(r(u),u)‖L22=1. By explicit computation, we get   r(u)=−1+1−‖u‖L22 (2.5) from which we deduce that $$r(u)$$ is well defined and smooth in a neighborhood of $$0$$ in $$H^1$$ and that moreover $$|r(u)|=\mathcal{O}{(\left\|{u}\right\|_{L^2}^2)}$$ when $$u$$ is small enough. Hence, $$u\mapsto\chi(r(u),u)$$ is a local parameterization of $$\mathcal{S} \cap V$$ in a neighborhood of $$\eta$$, where   S:={ψ∈V|‖ψ‖L2=1}. Note that in this parameterization, $$u = 0$$ corresponds to the ground state $$\eta$$. We then define   H(u):=H~(r(u),u)=H(χ(r(u),u)). (2.6) The main result in Weinstein (1985) (see also Grillakis et al., 1987, 1990 and Froehlich et al., 2004), is the following. Proposition 2.2 The point $$u=0$$ is a nondegenerate minimum of $${\mathcal{H}}$$: there exist some positive constants $$c_0$$ and $$\rho_0$$ such that   ∀v∈W,d2H(0).(v,v)≥c0‖v‖H12. (2.7) Note that we have   ∇ur(u)=−u1−‖u‖L22=−u+O(‖u‖H13). Hence, as $$ r(u) = \mathcal{O}(\left\|{u}\right\|_{H^1})$$, we have   ∇uH(u) =∂rH~(r(u),u)∇ur(u)+(∇uH~)(r(u),u) =λu−12Δu−3η2u+2⟨η3,u⟩η+O(‖u‖H12) =PW(λu−12Δu−3η2u)+O(‖u‖H12). From the previous proposition, we deduce the following. Corollary 2.3 The operator $$A: W \to W$$ defined by   Au:=PW(λu−12Δu−3η2u) (2.8) is $$L^2$$ symmetric and positive definite in $$H^1$$: we have   c‖u‖H12≤⟨u,Au⟩≤C‖u‖H12 (2.9) for some constants $$c$$ and $$C$$ and   ∇uH(u)=Au+O(‖u‖H12). (2.10) Let us remark that the coercitivity relation (2.9) combined with Cauchy–Schwarz inequality implies that   c‖u‖H12≤⟨u,Au⟩≤C‖u‖L2‖Au‖L2≤C‖u‖H1‖Au‖L2 for some constant $$C$$. We thus infer the existence of a constant $$c > 0$$ such that   ‖Au‖L2≥c‖u‖H1and‖Au‖L22≥c⟨u,Au⟩. (2.11) 3. Continuous normalized gradient flow We consider the continuous normalized gradient flow (see, for instance, Bao & Du, 2004)   ∂tψ=−∇H(ψ)+⟨∇H(ψ),ψ‖ψ‖L2⟩ψ‖ψ‖L2, (3.1) which is the projection of the standard gradient flow $$\partial_t\psi=-\nabla H(\psi)$$ onto the unit $$L^2$$ sphere. The local existence of an $$H^1$$ solution to (3.1) is guaranteed by a standard argument, using the fact that $${\it{\Delta}}$$ defines a semigroup in $$H^1$$, and that $$H^1$$ is an algebra in dimension $$1$$ (see (2.4)). 3.1 The gradient flow in local variables We assume for the moment that $$\psi(t)$$ remains in a ball sufficiently close to $$\eta$$ so that we can write $$\psi(t) = ( 1 + r(t)) \eta + u(t)$$, with $$r(t) > -1$$ and $$u(t) \in W$$. We have in this case $$P_\eta \psi = ( 1 + r) \eta$$ and $$P_W \psi = u$$. Note also that   ⟨∇H(ψ),ψ⟩ = ⟨Pη∇H,Pηψ⟩+⟨PW∇H,PWψ⟩ = (1+r)∂rH~+⟨∇uH~,u⟩ and that   ‖ψ‖L22=(1+r)2+‖u‖L22. It turns out that the relation $$\left\|{\psi}\right\|_{L^2}^2 = 1$$ implies that $$\left\|{u}\right\|_{L^2} \leq 1$$. In the following we will work only with functions $$u$$ satisfying this condition. Applying successively the operators $$P_\eta$$ and $$P_W$$ to equation (3.1), we obtain the relations   ∂tr = −∂rH~(r,u)+(1+r)2(1+r)2+‖u‖L22∂rH~+(1+r)(1+r)2+‖u‖L22⟨∇uH~,u⟩,∂tu = −∇uH~(r,u)+(1+r)u(1+r)2+‖u‖L22∂rH~+u(1+r)2+‖u‖L22⟨∇uH~,u⟩. As the $$L^2$$ norm of $$\psi$$ is preserved, we calculate directly that $$(1 + r)^2 + \left\|{u}\right\|_{L^2}^2$$ is preserved along the flow of this system (and is constant, equal to 1) and that we can solve $$r$$ in terms of $$u$$ as above (see (2.5)). We thus obtain the closed equation   ∂tu=−∇uH~(r(u),u)+(1+r(u))u∂rH~(r(u),u)+u⟨∇uH~,u⟩. (3.2) But we have   (∇uH~)(r(u),u)=∇uH(u)+∂rH~(r(u),u)u1−‖u‖L22; hence we can write equation (3.2) as   ∂tu = −∇uH(u)+u⟨u,∇uH(u)⟩−∂rH~(r(u),u)u1−‖u‖L22 +u1−‖u‖L22∂rH~(r(u),u)+∂rH~(r(u),u)u‖u‖L221−‖u‖L22, and hence   ∂tu=−∇uH(u)+u⟨u,∇uH(u)⟩, (3.3) which is the gradient flow in local coordinates $$(r(u),u)$$ on $$\mathcal{S}$$. Note that with the notation of the previous paragraph and using (2.10), we can write this equation as   ∂tu=−Au+R(u) (3.4) for some function $$R(u) = \mathcal{O}(\left\|{u}\right\|_{H^1}^2)$$ (see Definition 2.1). 3.2 Convergence of the flow We prove now that the solution of the normalized gradient flow (3.1) converges toward the ground state $$\eta$$ if the initial value is sufficiently close to it. Theorem 3.1 There exists $$\mu > 0$$ such that if $$u(t) \in W$$ is a solution of (3.3) with $$\left\|{u_0}\right\|_{H^1}$$ sufficiently small then we have   ‖u(t)‖H1≤C(u0)e−μt for all $$t >0$$ and for some constant $$C(u_0)$$ depending on $$u_0$$. Hence, if $$\psi(t) \in V \cap {\mathcal{S}}$$ is a solution of (3.1) such that $$\left\|{\psi(0)-\eta}\right\|_{H^1}$$ is small enough, then for all $$t$$,   ‖ψ(t)−η‖L2≤Ce−μt for some constant $$C$$ and $$\mu$$. Proof. As $$A$$ is symmetric, we calculate that   ∂t⟨u(t),Au(t)⟩=−2⟨Au(t),Au(t)⟩+2⟨R(u(t)),Au(t)⟩, where $$R(u)$$ is given by (3.4). Using (2.11), and the fact that for all $$u$$ and $$v$$ in $$H^1$$, we have $$\left\langle{u},{Av}\right\rangle \leq C \left\|{u}\right\|_{H^1} \left\|{v}\right\|_{H^1}$$, we obtain   |∂t⟨u(t),Au(t)⟩|≤−c⟨u(t),Au(t)⟩+2C‖R(u(t))‖H1‖u(t)‖H1 for some positive constants $$c$$ and $$C$$. As we have by definition $$\left\|{R(u(t))}\right\|_{H^1} \leq \left\|{u(t)}\right\|_{H^1}^2 \leq C \left\langle{u(t)},{A u(t)}\right\rangle$$, we infer the estimate   |∂t⟨u(t),Au(t)⟩|≤−c⟨u(t),Au(t)⟩+C⟨u(t),Au(t)⟩3/2. If $$\left\langle{u(t) },{Au(t)}\right\rangle \leq B^2$$ at $$ t= 0$$, we check that by a continuity argument, we have the rough estimate   ⟨u(t),Au(t)⟩≤B2e−(c−CB)t. The result easily follows by taking $$B$$ small enough and using (2.9). □ 4. Convergence result We consider the scheme described in Section 1, that from a function $$\psi_n$$ of $$L^2$$ norm $$1$$, close enough to $$\eta$$, we can write $$\psi_n = (1 + r(u_n)) \eta + u_n$$ with $$u_n \in W$$. We then define $$\psi_{n+1}$$ by the relations (1.5) and (1.6). Note that as $$\psi_{n+1}$$ is of $$L^2$$ norm $$1$$, we can define $$u_{n+1}$$ by the formula $$ \psi_{n+1} = (1 + r(u_{n+1})\eta + u_{n+1}$$ if $$\psi_{n+1}$$ is close enough to $${\it{\Gamma}}$$. The map $$u_n \mapsto u_{n+1}$$ satisfies the following relation. Lemma 4.1 Let $$R \leq \frac{1}{4}$$ be given. There exists $$\tau_0$$ such that for all $$\tau \leq \tau_0$$ the numerical scheme (1.5) and (1.6) is well defined from $$\mathcal{S}\cap \mathcal{U}(R)$$ to $$\mathcal{S}\cap \mathcal{U}(2R)$$ and is equivalent to an application $$u_n \mapsto u_{n+1}$$ satisfying   un+1=un+τPW(12Δun+1−λun+η2un+1+2η2un)+O(τ‖un,un+1‖H12), (4.1) which is a time discretization of the equation   ∂tu =PW(−λu+12Δu+3η2u)+O(‖u‖H12) = −Au+O(‖u‖H12), corresponding to the normalized gradient flow system (3.1). Proof. Note that with the assumption on $$R$$, we can assume that $$\left\|{u_n}\right\|_{H^1} \leq \tfrac{1}{2}$$. Let us calculate $$(r_n^*,u_n^*)$$ such that   ψn∗=χ(rn∗,un∗)=(1+rn∗)η+un∗ with $$u_n^* \in W$$, that is, $$\left\langle{\eta},{u_n^*}\right\rangle = 0$$. As $$r_n = r(u_n) = \mathcal{O}(\left\|{u_n}\right\|_{H^1}^2)$$, we have $$|\psi_n|^2 = \eta^2 + 2u_n \eta + \mathcal{O}(\left\|{u_n}\right\|_{H^1}^2)$$. Let us assume that $$|r_n^*| \leq 1$$; we can expand the Hamiltonian term   −∇H^(ψn,ψn∗) =12Δψn∗+|ψn|2ψn∗ =(1+rn∗)(−η3+λη)+12Δun∗+η2un∗+(η3+2unη2)(1+rn∗)+O(‖un,un∗‖H12), and hence we find   −∇H^(ψn,ψn∗)=λη(1+rn∗)+12Δun∗+η2un∗+2unη2(1+rn∗)+O(‖un,un∗‖H12). (4.2) By taking the projection $$P_\eta$$ of the equation (1.5), we obtain   1+rn∗=1+rn+τ(1+rn∗)λ+2τ⟨un,η3⟩(1+rn∗)+O(τ‖un,un∗‖H12). Remark 4.2 At this stage, the same calculations for the implicit–explicit or fully implicit schemes would yield a term depending on $$ \eta^3$$ in (4.2) which is neither in $$W$$ nor in $${\mathbb{R}} \eta$$. After taking a projection, this term would not vanish, while the ‘constant’ term $$\lambda \eta ( 1+ r_n^*)$$ is here orthogonal to $$W$$ (which reflects the fact that the ground state is exactly integrated by the linearly implicit scheme). We deduce that for $$\tau \leq \tau_0$$ small enough with respect to $$R$$ (and $$\lambda$$),   1+rn∗ =1+rn1−τλ−2τ⟨un,η3⟩+O(τ‖un,un∗‖H12) =11−τλ(1+rn+2τ1−τλ⟨un,η3⟩)+O(τ‖un,un∗‖H12), and by taking the projection on $$W$$, we obtain   un∗=un+τPW(12Δun∗+η2un∗+21−τλunη2)+O(τ‖un,un∗‖H12). (4.3) Let us now calculate an asymptotic expansion of the $$L^2$$ norm of $$\psi_n^*$$,   ‖ψn∗‖L22=(1+rn∗)2+‖un∗‖L22. We have   (1+rn∗)2 =(11−τλ)2(1+rn+2τ1−τλ⟨un,η3⟩)2+O(τ‖un,un∗‖H12) =(11−τλ)2((1+rn)2+4τ1−τλ⟨un,η3⟩)+O(τ‖un,un∗‖H12) and   ‖un∗‖L22=‖un‖L22+O(τ‖un,un∗‖H12). Hence, we calculate that   ‖ψn∗‖L2−1 =((1+rn∗)2+‖un∗‖L22)−1/2 =(1−τλ)((1+rn)2+4τ1−τλ⟨un,η3⟩+‖un‖L22+O(τ‖un,un∗‖H12))−1/2 =(1−τλ)(1+11−τλ4τ⟨un,η3⟩+O(τ‖un,un∗‖H12))−1/2 =1−τλ−2τ⟨un,η3⟩+O(τ‖un,un∗‖H12). Let us rewrite (4.3) as   un∗−τPW(12Δun∗+η2un∗)=un+τPW(21−τλunη2)+O(τ‖un,un∗‖H12). Hence, by multiplying by $$\left\|{\psi_n^*}\right\|_{L^2}^{-1}$$, we obtain   un+1−τPW(12Δun+1+η2un+1)=un−τλun+τPW(2unη2)+O(τ‖un,un+1‖H12), which shows formula (4.1). From this equation and using Definition 2.1, the implicit function theorem guarantees the local well-posedness of the application $$u_n \mapsto u_{n+1}$$ around $$0$$, with an estimate   ‖un+1‖H1≤2‖un‖H1 (4.4) uniformly with respect to $$\tau \leq \tau_0$$, with $$\tau_0$$ sufficiently small and provided $$u_n$$ is small enough in a bounded neighborhood of $$0$$. This yields the result. □ We can now prove the main result of this section. Theorem 4.3 There exist constants $$M$$, $$\tau_0$$, $$r$$ and $$C$$, such that if $$\psi_0 \in \mathcal{S} \cap \mathcal{U}(M)$$, then for all $$\tau \leq \tau_0$$ the solution of the numerical scheme (1.5) and (1.6) satisfies   ∀n,‖ψn−η‖H1≤Ce−rnτ. (4.5) Proof. Note that with the previous notation, if $$u_n$$ is small enough, we have $$\left\|{\psi_n - \eta}\right\|_{H^1} = \mathcal{O}(\left\|{u_n}\right\|_{H^1})$$. We will prove the following result, which implies (4.5) after a slight change of constants: there exist $$M$$, $$\tau_0$$, $$r$$ and $$C,$$ such that if $$\left\|{u_0}\right\|_{H^1} \leq M$$ and $$\tau \leq \tau_0$$ then every sequence $$(u_n)_{n \geq 0}$$ defined by (4.1) satisfies the estimate   ∀n,‖un‖H1≤Ce−rnτ. We rewrite system (4.1) as   un+1=un−τ(Lun+1+Bun)+O(τ‖un,un+1‖H12), with   Lu=PW(12Δu+η2u)andBu=PW(−λu+2η2u). We know that $$A = L + B$$ is coercive in $$H^1$$. Note moreover that $$B$$ is a symmetric bounded operator in $$W$$. We can write   un+1=un−τ(L+B)un+1+τB(un−un+1)+O(τ‖un,un+1‖H12). Now for $$\tau$$ sufficiently small, the operator $$I + \tau B$$ is symmetric, invertible, positive and bounded from $$W$$ to itself. Hence, we can rewrite the numerical scheme as   un+1=un−τ(I+τB)−1(L+B)un+1+O(τ‖un,un+1‖H12). Setting $$v_n = ( I + \tau B)^{1/2} u_n$$, we obtain the relation   vn+1=vn−τ(I+τB)−1/2(L+B)(I+τB)−1/2vn+1+O(τ‖vn,vn∗‖H12). In this new variable $$v_n$$, we can thus write the induction relation   vn+1=vn−τAτvn+1+τwn, (4.6) with $$ A_\tau = ( I + \tau B)^{-1/2}(L + B)( I + \tau B)^{-1/2}$$,   wn=O(‖vn,vn+1‖H12)andc‖v‖H12≤⟨Aτv,v⟩≤C‖v‖H12, (4.7) for some constants $$c$$ and $$C$$ independent of $$\tau$$. Note also that the inequality (2.11) holds for the operators $$A_\tau$$ uniformly in $$\tau \leq \tau_0$$ sufficiently small. From (4.6), we can write   ⟨Aτvn+τAτwn,vn+τwn⟩ =⟨Aτvn+1+τAτ2vn+1,vn+1+τAτvn+1⟩ =⟨Aτvn+1,vn+1⟩+2τ⟨Aτvn+1,Aτvn+1⟩+τ2⟨Aτ2vn+1,Aτvn+1⟩ ≥(1+2cτ)⟨Aτvn+1,vn+1⟩, using (2.11) and the fact that $$\left\langle{A_\tau^2 v_{n+1}},{A_\tau v_{n+1}}\right\rangle \geq 0$$. On the other hand, we have   ⟨Aτvn+τAτwn,vn+τwn⟩=⟨Aτvn,vn⟩+2τ⟨Aτvn,wn⟩+τ2⟨Aτwn,wn⟩. But using an integration by parts in the unbounded part of $$A_\tau$$, we have for all $$\tau \leq \tau_0$$ sufficiently small, and for any $$v,w$$ in $$H^1$$,   |⟨Aτv,w⟩|≤C‖v‖H1‖w‖H1. Hence we have, using the estimates (4.7) on $$w_n$$ (and Definition 2.1),   |⟨Aτvn,wn⟩|≤C⟨Aτvn,vn⟩3/2+C⟨Aτvn,vn⟩1/2⟨Aτvn+1,vn+1⟩ and   |⟨Aτwn,wn⟩|≤C⟨Aτvn,vn⟩2+C⟨Aτvn+1,vn+1⟩2. Gathering together the previous inequalities, we obtain that the sequence $$V_n = \left\langle{A_\tau v_n},{v_n}\right\rangle$$ satisfies the induction formula   (1+2cτ)Vn+1≤Vn+CτVn1/2Vn+1+CτVn3/2+Cτ2Vn2+Cτ2Vn+12 (4.8) and $$V_0 = M^2$$ (after a slight modification of the constant $$M$$). We will prove by induction that we can choose $$a$$ such that for all $$n \in {\mathbb{N}}$$,   Vn=⟨Aτvn,vn⟩≤M2e−anτ. Let us assume that this inequality holds for a given $$n$$. Using (4.4), we have that $$V_{n+1} \leq (C M)^2$$ for some constant $$C$$ independent of $$M$$ and $$\tau \leq \tau_0$$. Hence equation (4.8) implies that   (1+2cτ)Vn+1 ≤Vn+CτMVn+1+CτMVn+Cτ2M2Vn+Cτ2M2Vn+1 ≤ (1+CτM)Vn+CτMVn+1 after possible modifications of the constant $$C$$, and provided $$M$$ and $$\tau_0$$ are small enough. The previous inequality implies   Vn+1≤1+CMτ1+cτ−CMτVn+1. If $$M \tau_0$$ is small enough, we have for all $$\tau \leq \tau_0$$,   1+CMτ1+2cτ−CMτ≤e−aτ for some constant $$a$$ independent of $$\tau$$. This implies that $$V_{n+1} \leq {\it \,{e}}^{-a} V_n \leq M^2 e^{- a (n+1)}$$ and this concludes the induction argument. □ Remark 4.4 If we choose another discretization—implicit–explicit or fully implicit; see (1.7)—then as explained in Remark 4.2 the term in $$\eta^3$$ does not disappear in the equation, and after some manipulation we end up with a recursion formula of the form   vn+1=vn+τ2g−τAτvn+1+τρn with $$\rho_n = \mathcal{O}( \left\|{v_n,v_{n+1}}\right\|_{H^1}^2 ) $$ and $$g \in W$$ a nonzero function. It is relatively easy to prove that in this case, $$v_n$$ converges toward a function $$v_\tau^\infty = \tau A_{\tau}^{-1} g + \mathcal{O}(\tau^2)$$. Indeed, the previous equation is a discretization of the modified problem   ∂tv=τg−Aτv+B(v) whose solution exponentially converges toward the solution of the problem   τg−Aτv+B(v)=0, which exists by the standard implicit function theorem in $$H^1$$. This shows that in this case, $$\psi_n$$ converges exponentially in $$H^1$$ toward a modified soliton $$\eta_{\tau} = \eta + \mathcal{O}(\tau)$$. As the result is weaker than for the linearly implicit scheme, we do not give more mathematical details. 5. Fully discrete case The goal of this last section is to show that the previous results carry over to fully discrete systems. This fact relies essentially on the results in Bambusi & Penati (2010) and Bambusi et al. (2013). After discretization of the previous system by the finite difference method, there exists a discrete ground state $$\eta_h$$ minimizing a discrete convex functional $${\mathcal{H}}_h$$, which is an approximation of the exact functional $${\mathcal{H}}$$ defined above. Here, $$h$$ denotes the space discretization parameter and $$\eta_h$$ is close to $$\eta$$ with an error of order $$h$$. Moreover, the same holds true after Dirichlet cut-off upon adding an exponentially decreasing error term, and it can be proven that the discrete functionals satisfy the same estimates as the continuous ones, uniformly with respect to the discretization parameters. Having fixed a positive parameter $$h$$ we discretize in space by substituting the sequence $$\psi^{{\ell}}\simeq\psi(h {\ell})$$, $${\ell}\in{\mathbb{Z}}$$ for the function $$\psi(x)$$ and the second-order finite difference operator $${\it{\Delta}}_h$$ defined by   (Δhψ)ℓ:=ψℓ+1+ψℓ−1−2ψℓh2 (5.1) for the Laplace operator $$\partial_{xx}$$. We also impose Dirichlet boundary conditions for $$|j| \geq K+1$$, and the parabolic equation (1.4) then becomes   {ddtψℓ=12h2(ψℓ+1+ψℓ−1−2ψℓ)+|ψℓ|2ψℓ,ℓ∈Z,|ℓ|≤K,ψℓ=0,|ℓ|≥K+1,  (5.2) where $$t \mapsto \psi(t) = (\psi^\ell(t))_{\ell \in {\mathbb{Z}}}$$ is an application from $${\mathbb{R}}$$ to $${\mathbb{R}}^{\mathbb{Z}}$$. With this equation is associated a Hamiltonian function and a discrete $$L^2$$ norm given by   Hh(ψ)=h∑j∈Z[|ψj−ψj−1h|2−|ψj|42]andNh(ψ)=h∑j∈Z|ψj|2. (5.3) The discrete space of functions is   Vh,K={ψj∈CZ|ψj=ψ−j,ψj=0forj>K}, equipped with the discrete $$H^1$$ norm   ‖ψ‖h2=2h∑j∈Z|ψj+1−ψj|2h2+h∑j∈Z|ψj|2, (5.4) and we also define the scalar product   ⟨ψ,φ⟩h:=h∑j∈Zψjφj. Following Bambusi & Penati (2010) and Bambusi et al. (2013), we identify $$V_h$$ with a finite element subspace of $$H^1({\mathbb{R}})$$. More precisely, defining the function $$s:{\mathbb{R}} \to {\mathbb{R}}$$ by   s(x)={0if|x|>1,x+1if−1≤x≤0,−x+1if0≤x≤1,  (5.5) the identification is done through the map $$i_h: V_h \to H^1({\mathbb{R}})$$ defined by   {ψj}j∈Z↦(ihψ)(x):=∑j∈Zψjs(xh−j). (5.6) The main results in Bambusi & Penati (2010) and Bambusi et al. (2013) can expressed as follows. Theorem 5.1 For $$h$$ sufficiently small and $$K$$ sufficiently large, there exists a unique real minimizer $$\eta_{h,K} = (\eta_{h,k}^j)_{j \in {\mathbb{Z}}} \in V_{h,K}$$ of the functional $$H_h(\psi)$$ over the set   {ψ=(ψj)j∈Z∈Vh,K,Nh(ψ)=1} equipped with the norm (5.4). Moreover, $$\eta_{h,K}$$ satisfies the equation   12(Δhηh,K)ℓ+|ηh,Kℓ|2ηh,Kℓ =λhηh,Kℓ,ℓ=−K,…,K,ηh,Kℓ =0,|ℓ|≥K+1 (5.7) for some $$\lambda_h > 0$$. Moreover, we can define a local coordinate system $$\psi = \chi_h(r,u) = (1 + r) \eta_{h,K} + u$$ in a vicinity of $$\eta_{h,K}$$ in $$V_{h,K}$$ with $$u \in W_h := \{\hat{\mathrm{E}}\, u \in V_{h,K},\quad \left\langle{u},{\eta_{h,K}}\right\rangle_h = 0\}$$, and such that if $$r_h(u)$$ is defined by the implicit relation $$N_h( \chi_h(r_h(u),u))) = 1$$ then the functional $${\mathcal{H}}_{h}(u) = H_h(\chi_h(r_h(u),u)))$$ has a unique nondegenerate minimum in $$u = 0$$ in $$V_{h,K}$$ and satisfies a convexity estimate of the form (2.7) uniformly in $$h$$ and $$K$$. Finally, the discrete ground state $$\eta_{h,K}$$ is close to the continuous one in the sense that   ‖ihηh,K−η‖H1≤C(h+1h2e−C1hK), (5.8) where $$C$$ and $$C_1$$ do not depend on $$h$$ and where $$i_h$$ is the embedding (5.6). The fully discrete algorithm corresponding to (1.5) and (1.6) then consists of the following two steps. From $$(\psi^\ell_n)_{j = -K}^K$$, such that $$N_h(\psi_n) = 1$$, (i) compute $$\psi^*_n = (\psi^{*,\ell}_n)_{j \in {\mathbb{Z}}}$$ defined by the relation   {ψn∗,ℓ=ψnℓ+τ(12(Δhψn∗)ℓ+|ψnℓ|2ψn∗,ℓ),ℓ∈Z,|ℓ|≤K,ψ∗,ℓ=0,|ℓ|≥K+1;  (5.9) (ii) normalize:   ψn+1ℓ=ψn∗,ℓNh(ψn∗). (5.10) Equation (5.7) guarantees that this fully discrete algorithm preserves exactly the discrete ground state $$\eta_{h,K}$$. We can now copy the proof of the continuous case and adapt it directly to the fully discrete case: we can prove that   ‖ψn−ηh,K‖h≤Ce−cnτ,n≥0, for some constants $$C$$ and $$c$$ independent of $$h$$, $$K$$ and $$\tau$$. Using the previous estimate, we obtain the following result. Theorem 5.2 There exist constants $$M$$, $$C$$, $$C_1$$ and $$c$$ such that for $$h$$, $$\tau$$ sufficiently small and $$K$$ sufficiently large, if $$(\psi_0) = (\psi^\ell_0)_{\ell = -K}^{K}$$ satisfies   ‖ihψ0−η‖H1≤M then the solution $$\psi_n$$ of the fully discrete algorithm (5.9) and (5.10) satisfies   ‖ihψn−η‖H1≤C(e−cnτ+h+1h2e−C1hK). Remark 5.3 Following Remark 4.4, we can prove a similar result for the implicit–explicit and fully implicit algorithms (1.7), but the fully discrete schemes will converge toward a discrete ground state $$\eta_{\tau,h,K}$$ satisfying   ‖ihητ,h,K−η‖H1≤C(τ+h+1h2e−C1hK), and the previous result has to be modified accordingly in these cases. Funding ERC Starting Grant GEOPARDI (No. 279389). References Adhikari S. K. ( 2000a) Numerical solution of the two-dimensional Gross-Pitaevskii equation for trapped interacting atoms. Phys. Lett. A , 265, 91– 96. Google Scholar CrossRef Search ADS   Adhikari S. K. ( 2000b) Numerical solution of spherically symmetric Gross-Pitaevskii equation in two space dimensions. Phys. Lett. E , 62, 2937– 2944. Bambusi D., Faou E. & Grébert B. ( 2013) Existence and stability of ground states for fully discrete approximations of the nonlinear Schrödinger equation. Numer. Math. , 123, 461– 492. Google Scholar CrossRef Search ADS   Bambusi D. & Penati T. ( 2010) Continuous approximation of breathers in one and two dimensional DNLS lattices. Nonlinearity , 23, 143– 157. Google Scholar CrossRef Search ADS   Bao W. & Du Q. ( 2004) Computing the ground state solution of Bose-Einstein condensates by a normalized gradient flow. SIAM J. Sci. Comput. , 25, 1674– 1697. Google Scholar CrossRef Search ADS   Bao W. & Tang W. ( 2003) Ground state solution of trapped interacting Bose-Einstein condensate by directly minimizing the energy functional. J. Comput. Phys. , 187, 230– 254. Google Scholar CrossRef Search ADS   Chiofalo M. L., Succi S. & Tosi M. P. ( 2000) Ground state of trapped interacting Bose-Einstein condensates by an explicit imaginary time algorithm. Phys. Rev. E , 62, 7438– 7444. Google Scholar CrossRef Search ADS   Edwards M. & Burnett K. ( 1995) Numerical solution of the nonlinear Schrödinger equation for small samples of trapped neutral atoms. Phys. Rev. , 1, 1382– 1386. Google Scholar CrossRef Search ADS   Froehlich J., Gustafson S., Jonsson L. & Sigal I. M. ( 2004) Solitary wave dynamics in an external potential. Comm. Math. Phys. , 250, 613– 642. Google Scholar CrossRef Search ADS   Grillakis M., Shatah H. & Strauss W. ( 1987) Stability theory of solitary waves in the presence of symmetry. I. J. Funct. Anal. , 74, 160– 197. Google Scholar CrossRef Search ADS   Grillakis M., Shatah H. & Strauss W. ( 1990) Stability theory of solitary waves in the presence of symmetry. II. J. Funct. Anal. , 94, 308– 348. Google Scholar CrossRef Search ADS   Hairer E. & Lubich C. ( 2014) Energy-diminishing integration of gradient systems. IMA J. Numer. Anal. , 34, 452– 461. Google Scholar CrossRef Search ADS   Weinstein M. I. ( 1985) Modulational stability of ground states of nonlinear Schrödinger equations. SIAM J. Math. Anal. , 16, 472– 491. Google Scholar CrossRef Search ADS   © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Numerical Analysis Oxford University Press

Convergence of a normalized gradient algorithm for computing ground states

Download PDF
17 pages

Loading next page...
 
/lp/ou_press/convergence-of-a-normalized-gradient-algorithm-for-computing-ground-jDw6DKURDE

References (14)

Publisher
Oxford University Press
Copyright
© The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.
ISSN
0272-4979
eISSN
1464-3642
DOI
10.1093/imanum/drx009
Publisher site
See Article on Publisher Site

Abstract

Abstract We consider the approximation of the ground state of the one-dimensional cubic nonlinear Schrödinger equation by a normalized gradient algorithm combined with a linearly implicit time integrator, and a finite difference space approximation. We show that this method, also called the imaginary time evolution method in the physics literature, is locally convergent, and we provide error estimates: for initial data in a neighborhood of the ground state, the algorithm converges exponentially toward a modified soliton that is a space discretization of the exact soliton, with error estimates depending on the discretization parameters. 1. Introduction The goal of this article is to give a convergence proof for a normalized gradient algorithm used to numerically compute ground states of Schrödinger equations fulfilling symmetry and coercivity conditions as considered in the seminal works of Weinstein (1985) and Grillakis et al. (1987, 1990). This algorithm is also called the imaginary time method in the physics literature; see for instance Edwards & Burnett (1995), Chiofalo et al. (2000), Adhikari (2000a,b), Bao & Tang (2003) and the references therein. Let us describe the algorithm in the case of the focusing cubic nonlinear Schrödinger equation   i∂tψ=−12Δψ−|ψ|2ψ, (NLS) set on $${\mathbb{R}}$$, where $$\psi(t,x)$$ depends on space variables $$x \in {\mathbb{R}}$$. With this equation is associated the energy   H(ψ,ψ¯)=14∫R|∇ψ|2−|ψ|4, (1.1) which is preserved by the flow of (NLS) for all times. The equation (NLS) can be written as   i∂tψ=−12Δψ−|ψ|2ψ=2∂H∂ψ¯(ψ,ψ¯). In the rest of this article, the notation $$\nabla H$$ will denote the $$L^2$$ derivative of the energy $$H$$ with respect to real functions $$\psi$$. Note that we have for a real function $$u$$,   ∇H(u)=2∂H∂ψ¯(u,u)=−12Δu−u3, the left-hand side denoting the Fréchet derivative of $$H(u)$$ considered as a functional acting on real functions. Note that naturally, $$\nabla H \in H^{-1}$$ with the embedding $$H^{-1} \subset L^2 \subset H^1$$. With this notation, the ground state $$\eta(x)$$ is defined as the unique real symmetric minimizer in $$H^1$$ of the problem   min‖ψ‖L2=1H(ψ). (1.2) In the one-dimensional cubic case considered in this article, explicit computations show that   η(x):=12sech(x2). We denote by $$\lambda$$ the Lagrange multiplier associated with this minimization problem, such that   ∇H(η)=−12Δη−η3=−λη. (1.3) With the ground state $$\eta$$ is associated the solution $$\eta(t,x) = \eta(x) e^{{\rm i} \lambda t}$$ of the time-dependent equation (NLS). Using translation and scaling, the ground state gives rise to a family of explicit solutions of (NLS) that have the property of being orbitally stable in $$H^1$$ (see Weinstein, 1985; Grillakis et al., 1987, 1990; Bambusi et al., 2013). For instance, any solution starting close to $$\eta(x)$$ will remain close to the manifold $$\{\,{\it e}^{{\rm i}\alpha} \eta ( x - c)\, | \alpha,c \in {\mathbb{R}}\} \subset H^1$$, the rotation and translation being natural invariant group actions of the nonlinear Schrödinger equation (NLS). In more general situations, $$\eta$$ is not explicitly known, and one has to rely on numerical simulations to compute it. To this aim, the imaginary time method, which is a nonlinear version of the normalized gradient algorithm, is widely used. The goal of the present article is to analyse the efficiency of this method in the simple case described above (the results obtained are in fact valid in more general situations—essentially all situations where the Grillakis–Shatah–Strauss arguments apply). The time-discretized algorithm consists in defining a sequence a functions $$\{\psi_n\}_{n \in {\mathbb{N}}}$$ as follows: (i) For a given $$\psi_n$$, an intermediate function $$\psi_{n}^*$$ is defined as a numerical approximation of the solution of the parabolic equation   ∂tψ=12Δψ+|ψ|2ψ=−∇H(ψ), (1.4) starting in $$\psi_n$$, over a time interval $$[0,\tau]$$, where $$\tau$$ is a given time step. To compute $$\psi_n^*$$, we will see that the best results are given by the linearly implicit method   ψn∗=ψn−τ∇H^(ψn,ψn∗), (1.5) where   −∇H^(ψn,ψn∗)=12Δψn∗+|ψn|2ψn∗. (ii) Then, we define the normalized function   ψn+1=ψn∗‖ψn∗‖L2. (1.6) Note that the algorithm presented above preserves the real nature of the function $$\psi$$: if $$\psi_0$$ is real valued then for all $$n \in {\mathbb{N}}$$, $$\psi_n$$ is real, and the same holds true for symmetric functions satisfying $$\psi_n(-x) = \psi_n(x)$$. Our results can be summarized as follows: In the semidiscrete case described above, if the initial data $$\psi_0$$ is real, symmetric and sufficiently close to $$\eta$$ then $$\psi_n$$ converges exponentially toward $$\eta$$ in $$H^1$$. In the fully discrete case, where the discretization in space is made by finite differences in space with mesh $$h$$ combined with a Dirichlet cut-off for large values of $$x$$ (namely $$x \leq Kh$$), then we prove exponential convergence toward a modified soliton $$\eta_{h,K}$$ that is $$\mathcal{O}( h + \frac{1}{h^2} {\it \,{e}}^{- C_1 Kh})$$ close to the exact soliton $$\eta$$. The main property explaining the excellent performance of this method is the fact that the linearly implicit numerical scheme exactly preserves the ground state: if $$\psi_n = \eta$$ then $$\psi_{n+1} = \eta$$, and the same holds true for discrete-in-space ground states satisfying a discrete version of (1.3). This fact is very general: it holds for any linearly implicit scheme applied to a semilinear partial differential equation (PDE). This important feature of course makes this scheme much more attractive than other possible schemes, where the nonlinearity would be approached by   −∇H^(ψn,ψn∗)=12Δψn∗+{|ψn|2ψn(semi−explicit),|ψn∗|2ψn∗(fullyimplicit),  (1.7) despite the fact that fully implicit schemes enjoy the energy-diminishing property for a standard gradient system (see Hairer & Lubich, 2014). As we will discuss later, these schemes still converge, but toward modified ground states $$\eta_\tau = \eta + \mathcal{O}(\tau)$$. As these results would be weaker, with longer detailed proofs, we give only the main arguments to obtain them, keeping a full detailed convergence proof for the linearly implicit scheme. 2. The Hamiltonian near the ground state We work in the space $$V$$ of real symmetric functions of $$H^1$$:   V:={φ∈H1(R,R)|φ(−x)=φ(x)}. We consider the usual $$L^2$$ and $$H^1$$ norms and denote by $$\left\langle{\cdot},{\cdot}\right\rangle$$ the canonical real scalar product on $$L^2$$:   ⟨φ,ψ⟩=∫Rφ(x)ψ(x)ds,‖φ‖L22=⟨φ,φ⟩,‖φ‖H12=‖φ‖L22+‖∂xφ‖L22. 2.1 Coordinates in the neighborhood of the ground state We introduce the set of the functions $$R$$-close to $$\eta$$,   U(R):={φ∈V|‖φ−η‖H1<R}, the set   W:={u∈V|⟨u,η⟩=0} and the map $$\chi$$  χ:R×W→V,(r,u)↦(1+r)η+u. The map $$\chi$$ allows us to use $$(r,u)$$ as coordinates in $$\,\mathcal{U}(R)$$. Observe that $$\chi$$ is smooth with bounded derivatives and so is the inverse $$\chi^{-1}$$, from the explicit formula   χ−1:V→R×W,ψ↦(r(ψ),u(ψ))=(⟨ψ,η⟩−1,ψ−⟨ψ,η⟩η). We will also use the following notation: we define the $$L^2$$ projectors   Pηu=⟨u,η⟩ηandPWu=I−Pη, and we define the function   H~(r,u)=(H∘χ)(r,u). (2.1) We can verify the following relations:   ∂rH~(r,u)=⟨∇H(χ(r,u)),η⟩=η−1(Pη∇H)(χ(r,u))∈R (2.2) and   ∇uH~(r,u)=PW∇H(χ(r,u))∈H−1. (2.3) Before collecting some expressions of the Hamiltonian and the gradient flow in coordinates $$(r,u)$$, let us introduce the following notation. Definition 2.1 Let $$p \geq 1$$. We say that $$R(u) = \mathcal{O}(\left\|{u}\right\|_{H^1}^p)$$ if $$R: H^1 \to H^1$$ is a smooth function such that for all $$B$$, there exists a constant $$C$$ such that for all $$u\in H^1$$, $$\left\|{u}\right\|_{H^1} \leq B$$, we have   ‖R(u)‖H1≤C‖u‖H1p. We will also use the notation $$R(u,v) = \mathcal{O}(\left\|{u,v}\right\|_{H^1}^p)$$ for a smooth function $$R: H^1 \times H^1 \to H^1$$ if for all $$B$$, there exists $$C$$ such that for all $$u,v \in H^1$$ satisfying $$\left\|{u}\right\|_{H^1} \leq B$$ and $$\left\|{v}\right\|_{H^1} \leq B$$ then we have   ‖R(u,v)‖H1≤C(‖u‖H1p+‖v‖H1p). Finally, a function $$u = \mathcal{O}(r)$$ if $$\left\|{u}\right\|_{H^1} \leq C r$$ for $$r$$ small enough and a constant $$C$$ independent of $$u$$. With this notation and the fact that $$\psi = (1 + r) \eta + u$$, we compute using (1.3) that   −∇H(ψ) =12Δψ+ψ3 = (1+r)(−η3+λη)+12Δu+u3 +(1+r)3η3+3(1+r)2η2u+3(1+r)ηu2 =λη+12Δu+3η2u+O(r+‖u‖H12). Note that to obtain the bound, we have used the fact that $$\eta \in H^1$$, as well as the estimate   ‖uv‖H1≤C‖u‖H1‖v‖H1 (2.4) for two functions $$u$$ and $$v$$. We deduce using (2.2) and (2.3) that   ∂rH~(r,u) = −λ−12⟨u,Δη⟩−3⟨η3,u⟩+O(r+‖u‖H12) = −λ−2⟨η3,u⟩+O(r+‖u‖H12), as $$\left\langle{u},{\eta}\right\rangle = 0$$ and $$\left\langle{\eta},{\eta}\right\rangle = 1$$ and   ∇uH~(r,u) = −PW(λη+12Δu+3η2u)+O(r+‖u‖H12) = −12Δu−3η2u+2⟨η3,u⟩η+O(r+‖u‖H12). 2.2 Projection onto the unit $$L^2$$ sphere Let us now define the function $$u\mapsto r(u)$$ from $$W$$ to $${\mathbb{R}}$$ by the implicit relation   ‖χ(r(u),u)‖L22=1. By explicit computation, we get   r(u)=−1+1−‖u‖L22 (2.5) from which we deduce that $$r(u)$$ is well defined and smooth in a neighborhood of $$0$$ in $$H^1$$ and that moreover $$|r(u)|=\mathcal{O}{(\left\|{u}\right\|_{L^2}^2)}$$ when $$u$$ is small enough. Hence, $$u\mapsto\chi(r(u),u)$$ is a local parameterization of $$\mathcal{S} \cap V$$ in a neighborhood of $$\eta$$, where   S:={ψ∈V|‖ψ‖L2=1}. Note that in this parameterization, $$u = 0$$ corresponds to the ground state $$\eta$$. We then define   H(u):=H~(r(u),u)=H(χ(r(u),u)). (2.6) The main result in Weinstein (1985) (see also Grillakis et al., 1987, 1990 and Froehlich et al., 2004), is the following. Proposition 2.2 The point $$u=0$$ is a nondegenerate minimum of $${\mathcal{H}}$$: there exist some positive constants $$c_0$$ and $$\rho_0$$ such that   ∀v∈W,d2H(0).(v,v)≥c0‖v‖H12. (2.7) Note that we have   ∇ur(u)=−u1−‖u‖L22=−u+O(‖u‖H13). Hence, as $$ r(u) = \mathcal{O}(\left\|{u}\right\|_{H^1})$$, we have   ∇uH(u) =∂rH~(r(u),u)∇ur(u)+(∇uH~)(r(u),u) =λu−12Δu−3η2u+2⟨η3,u⟩η+O(‖u‖H12) =PW(λu−12Δu−3η2u)+O(‖u‖H12). From the previous proposition, we deduce the following. Corollary 2.3 The operator $$A: W \to W$$ defined by   Au:=PW(λu−12Δu−3η2u) (2.8) is $$L^2$$ symmetric and positive definite in $$H^1$$: we have   c‖u‖H12≤⟨u,Au⟩≤C‖u‖H12 (2.9) for some constants $$c$$ and $$C$$ and   ∇uH(u)=Au+O(‖u‖H12). (2.10) Let us remark that the coercitivity relation (2.9) combined with Cauchy–Schwarz inequality implies that   c‖u‖H12≤⟨u,Au⟩≤C‖u‖L2‖Au‖L2≤C‖u‖H1‖Au‖L2 for some constant $$C$$. We thus infer the existence of a constant $$c > 0$$ such that   ‖Au‖L2≥c‖u‖H1and‖Au‖L22≥c⟨u,Au⟩. (2.11) 3. Continuous normalized gradient flow We consider the continuous normalized gradient flow (see, for instance, Bao & Du, 2004)   ∂tψ=−∇H(ψ)+⟨∇H(ψ),ψ‖ψ‖L2⟩ψ‖ψ‖L2, (3.1) which is the projection of the standard gradient flow $$\partial_t\psi=-\nabla H(\psi)$$ onto the unit $$L^2$$ sphere. The local existence of an $$H^1$$ solution to (3.1) is guaranteed by a standard argument, using the fact that $${\it{\Delta}}$$ defines a semigroup in $$H^1$$, and that $$H^1$$ is an algebra in dimension $$1$$ (see (2.4)). 3.1 The gradient flow in local variables We assume for the moment that $$\psi(t)$$ remains in a ball sufficiently close to $$\eta$$ so that we can write $$\psi(t) = ( 1 + r(t)) \eta + u(t)$$, with $$r(t) > -1$$ and $$u(t) \in W$$. We have in this case $$P_\eta \psi = ( 1 + r) \eta$$ and $$P_W \psi = u$$. Note also that   ⟨∇H(ψ),ψ⟩ = ⟨Pη∇H,Pηψ⟩+⟨PW∇H,PWψ⟩ = (1+r)∂rH~+⟨∇uH~,u⟩ and that   ‖ψ‖L22=(1+r)2+‖u‖L22. It turns out that the relation $$\left\|{\psi}\right\|_{L^2}^2 = 1$$ implies that $$\left\|{u}\right\|_{L^2} \leq 1$$. In the following we will work only with functions $$u$$ satisfying this condition. Applying successively the operators $$P_\eta$$ and $$P_W$$ to equation (3.1), we obtain the relations   ∂tr = −∂rH~(r,u)+(1+r)2(1+r)2+‖u‖L22∂rH~+(1+r)(1+r)2+‖u‖L22⟨∇uH~,u⟩,∂tu = −∇uH~(r,u)+(1+r)u(1+r)2+‖u‖L22∂rH~+u(1+r)2+‖u‖L22⟨∇uH~,u⟩. As the $$L^2$$ norm of $$\psi$$ is preserved, we calculate directly that $$(1 + r)^2 + \left\|{u}\right\|_{L^2}^2$$ is preserved along the flow of this system (and is constant, equal to 1) and that we can solve $$r$$ in terms of $$u$$ as above (see (2.5)). We thus obtain the closed equation   ∂tu=−∇uH~(r(u),u)+(1+r(u))u∂rH~(r(u),u)+u⟨∇uH~,u⟩. (3.2) But we have   (∇uH~)(r(u),u)=∇uH(u)+∂rH~(r(u),u)u1−‖u‖L22; hence we can write equation (3.2) as   ∂tu = −∇uH(u)+u⟨u,∇uH(u)⟩−∂rH~(r(u),u)u1−‖u‖L22 +u1−‖u‖L22∂rH~(r(u),u)+∂rH~(r(u),u)u‖u‖L221−‖u‖L22, and hence   ∂tu=−∇uH(u)+u⟨u,∇uH(u)⟩, (3.3) which is the gradient flow in local coordinates $$(r(u),u)$$ on $$\mathcal{S}$$. Note that with the notation of the previous paragraph and using (2.10), we can write this equation as   ∂tu=−Au+R(u) (3.4) for some function $$R(u) = \mathcal{O}(\left\|{u}\right\|_{H^1}^2)$$ (see Definition 2.1). 3.2 Convergence of the flow We prove now that the solution of the normalized gradient flow (3.1) converges toward the ground state $$\eta$$ if the initial value is sufficiently close to it. Theorem 3.1 There exists $$\mu > 0$$ such that if $$u(t) \in W$$ is a solution of (3.3) with $$\left\|{u_0}\right\|_{H^1}$$ sufficiently small then we have   ‖u(t)‖H1≤C(u0)e−μt for all $$t >0$$ and for some constant $$C(u_0)$$ depending on $$u_0$$. Hence, if $$\psi(t) \in V \cap {\mathcal{S}}$$ is a solution of (3.1) such that $$\left\|{\psi(0)-\eta}\right\|_{H^1}$$ is small enough, then for all $$t$$,   ‖ψ(t)−η‖L2≤Ce−μt for some constant $$C$$ and $$\mu$$. Proof. As $$A$$ is symmetric, we calculate that   ∂t⟨u(t),Au(t)⟩=−2⟨Au(t),Au(t)⟩+2⟨R(u(t)),Au(t)⟩, where $$R(u)$$ is given by (3.4). Using (2.11), and the fact that for all $$u$$ and $$v$$ in $$H^1$$, we have $$\left\langle{u},{Av}\right\rangle \leq C \left\|{u}\right\|_{H^1} \left\|{v}\right\|_{H^1}$$, we obtain   |∂t⟨u(t),Au(t)⟩|≤−c⟨u(t),Au(t)⟩+2C‖R(u(t))‖H1‖u(t)‖H1 for some positive constants $$c$$ and $$C$$. As we have by definition $$\left\|{R(u(t))}\right\|_{H^1} \leq \left\|{u(t)}\right\|_{H^1}^2 \leq C \left\langle{u(t)},{A u(t)}\right\rangle$$, we infer the estimate   |∂t⟨u(t),Au(t)⟩|≤−c⟨u(t),Au(t)⟩+C⟨u(t),Au(t)⟩3/2. If $$\left\langle{u(t) },{Au(t)}\right\rangle \leq B^2$$ at $$ t= 0$$, we check that by a continuity argument, we have the rough estimate   ⟨u(t),Au(t)⟩≤B2e−(c−CB)t. The result easily follows by taking $$B$$ small enough and using (2.9). □ 4. Convergence result We consider the scheme described in Section 1, that from a function $$\psi_n$$ of $$L^2$$ norm $$1$$, close enough to $$\eta$$, we can write $$\psi_n = (1 + r(u_n)) \eta + u_n$$ with $$u_n \in W$$. We then define $$\psi_{n+1}$$ by the relations (1.5) and (1.6). Note that as $$\psi_{n+1}$$ is of $$L^2$$ norm $$1$$, we can define $$u_{n+1}$$ by the formula $$ \psi_{n+1} = (1 + r(u_{n+1})\eta + u_{n+1}$$ if $$\psi_{n+1}$$ is close enough to $${\it{\Gamma}}$$. The map $$u_n \mapsto u_{n+1}$$ satisfies the following relation. Lemma 4.1 Let $$R \leq \frac{1}{4}$$ be given. There exists $$\tau_0$$ such that for all $$\tau \leq \tau_0$$ the numerical scheme (1.5) and (1.6) is well defined from $$\mathcal{S}\cap \mathcal{U}(R)$$ to $$\mathcal{S}\cap \mathcal{U}(2R)$$ and is equivalent to an application $$u_n \mapsto u_{n+1}$$ satisfying   un+1=un+τPW(12Δun+1−λun+η2un+1+2η2un)+O(τ‖un,un+1‖H12), (4.1) which is a time discretization of the equation   ∂tu =PW(−λu+12Δu+3η2u)+O(‖u‖H12) = −Au+O(‖u‖H12), corresponding to the normalized gradient flow system (3.1). Proof. Note that with the assumption on $$R$$, we can assume that $$\left\|{u_n}\right\|_{H^1} \leq \tfrac{1}{2}$$. Let us calculate $$(r_n^*,u_n^*)$$ such that   ψn∗=χ(rn∗,un∗)=(1+rn∗)η+un∗ with $$u_n^* \in W$$, that is, $$\left\langle{\eta},{u_n^*}\right\rangle = 0$$. As $$r_n = r(u_n) = \mathcal{O}(\left\|{u_n}\right\|_{H^1}^2)$$, we have $$|\psi_n|^2 = \eta^2 + 2u_n \eta + \mathcal{O}(\left\|{u_n}\right\|_{H^1}^2)$$. Let us assume that $$|r_n^*| \leq 1$$; we can expand the Hamiltonian term   −∇H^(ψn,ψn∗) =12Δψn∗+|ψn|2ψn∗ =(1+rn∗)(−η3+λη)+12Δun∗+η2un∗+(η3+2unη2)(1+rn∗)+O(‖un,un∗‖H12), and hence we find   −∇H^(ψn,ψn∗)=λη(1+rn∗)+12Δun∗+η2un∗+2unη2(1+rn∗)+O(‖un,un∗‖H12). (4.2) By taking the projection $$P_\eta$$ of the equation (1.5), we obtain   1+rn∗=1+rn+τ(1+rn∗)λ+2τ⟨un,η3⟩(1+rn∗)+O(τ‖un,un∗‖H12). Remark 4.2 At this stage, the same calculations for the implicit–explicit or fully implicit schemes would yield a term depending on $$ \eta^3$$ in (4.2) which is neither in $$W$$ nor in $${\mathbb{R}} \eta$$. After taking a projection, this term would not vanish, while the ‘constant’ term $$\lambda \eta ( 1+ r_n^*)$$ is here orthogonal to $$W$$ (which reflects the fact that the ground state is exactly integrated by the linearly implicit scheme). We deduce that for $$\tau \leq \tau_0$$ small enough with respect to $$R$$ (and $$\lambda$$),   1+rn∗ =1+rn1−τλ−2τ⟨un,η3⟩+O(τ‖un,un∗‖H12) =11−τλ(1+rn+2τ1−τλ⟨un,η3⟩)+O(τ‖un,un∗‖H12), and by taking the projection on $$W$$, we obtain   un∗=un+τPW(12Δun∗+η2un∗+21−τλunη2)+O(τ‖un,un∗‖H12). (4.3) Let us now calculate an asymptotic expansion of the $$L^2$$ norm of $$\psi_n^*$$,   ‖ψn∗‖L22=(1+rn∗)2+‖un∗‖L22. We have   (1+rn∗)2 =(11−τλ)2(1+rn+2τ1−τλ⟨un,η3⟩)2+O(τ‖un,un∗‖H12) =(11−τλ)2((1+rn)2+4τ1−τλ⟨un,η3⟩)+O(τ‖un,un∗‖H12) and   ‖un∗‖L22=‖un‖L22+O(τ‖un,un∗‖H12). Hence, we calculate that   ‖ψn∗‖L2−1 =((1+rn∗)2+‖un∗‖L22)−1/2 =(1−τλ)((1+rn)2+4τ1−τλ⟨un,η3⟩+‖un‖L22+O(τ‖un,un∗‖H12))−1/2 =(1−τλ)(1+11−τλ4τ⟨un,η3⟩+O(τ‖un,un∗‖H12))−1/2 =1−τλ−2τ⟨un,η3⟩+O(τ‖un,un∗‖H12). Let us rewrite (4.3) as   un∗−τPW(12Δun∗+η2un∗)=un+τPW(21−τλunη2)+O(τ‖un,un∗‖H12). Hence, by multiplying by $$\left\|{\psi_n^*}\right\|_{L^2}^{-1}$$, we obtain   un+1−τPW(12Δun+1+η2un+1)=un−τλun+τPW(2unη2)+O(τ‖un,un+1‖H12), which shows formula (4.1). From this equation and using Definition 2.1, the implicit function theorem guarantees the local well-posedness of the application $$u_n \mapsto u_{n+1}$$ around $$0$$, with an estimate   ‖un+1‖H1≤2‖un‖H1 (4.4) uniformly with respect to $$\tau \leq \tau_0$$, with $$\tau_0$$ sufficiently small and provided $$u_n$$ is small enough in a bounded neighborhood of $$0$$. This yields the result. □ We can now prove the main result of this section. Theorem 4.3 There exist constants $$M$$, $$\tau_0$$, $$r$$ and $$C$$, such that if $$\psi_0 \in \mathcal{S} \cap \mathcal{U}(M)$$, then for all $$\tau \leq \tau_0$$ the solution of the numerical scheme (1.5) and (1.6) satisfies   ∀n,‖ψn−η‖H1≤Ce−rnτ. (4.5) Proof. Note that with the previous notation, if $$u_n$$ is small enough, we have $$\left\|{\psi_n - \eta}\right\|_{H^1} = \mathcal{O}(\left\|{u_n}\right\|_{H^1})$$. We will prove the following result, which implies (4.5) after a slight change of constants: there exist $$M$$, $$\tau_0$$, $$r$$ and $$C,$$ such that if $$\left\|{u_0}\right\|_{H^1} \leq M$$ and $$\tau \leq \tau_0$$ then every sequence $$(u_n)_{n \geq 0}$$ defined by (4.1) satisfies the estimate   ∀n,‖un‖H1≤Ce−rnτ. We rewrite system (4.1) as   un+1=un−τ(Lun+1+Bun)+O(τ‖un,un+1‖H12), with   Lu=PW(12Δu+η2u)andBu=PW(−λu+2η2u). We know that $$A = L + B$$ is coercive in $$H^1$$. Note moreover that $$B$$ is a symmetric bounded operator in $$W$$. We can write   un+1=un−τ(L+B)un+1+τB(un−un+1)+O(τ‖un,un+1‖H12). Now for $$\tau$$ sufficiently small, the operator $$I + \tau B$$ is symmetric, invertible, positive and bounded from $$W$$ to itself. Hence, we can rewrite the numerical scheme as   un+1=un−τ(I+τB)−1(L+B)un+1+O(τ‖un,un+1‖H12). Setting $$v_n = ( I + \tau B)^{1/2} u_n$$, we obtain the relation   vn+1=vn−τ(I+τB)−1/2(L+B)(I+τB)−1/2vn+1+O(τ‖vn,vn∗‖H12). In this new variable $$v_n$$, we can thus write the induction relation   vn+1=vn−τAτvn+1+τwn, (4.6) with $$ A_\tau = ( I + \tau B)^{-1/2}(L + B)( I + \tau B)^{-1/2}$$,   wn=O(‖vn,vn+1‖H12)andc‖v‖H12≤⟨Aτv,v⟩≤C‖v‖H12, (4.7) for some constants $$c$$ and $$C$$ independent of $$\tau$$. Note also that the inequality (2.11) holds for the operators $$A_\tau$$ uniformly in $$\tau \leq \tau_0$$ sufficiently small. From (4.6), we can write   ⟨Aτvn+τAτwn,vn+τwn⟩ =⟨Aτvn+1+τAτ2vn+1,vn+1+τAτvn+1⟩ =⟨Aτvn+1,vn+1⟩+2τ⟨Aτvn+1,Aτvn+1⟩+τ2⟨Aτ2vn+1,Aτvn+1⟩ ≥(1+2cτ)⟨Aτvn+1,vn+1⟩, using (2.11) and the fact that $$\left\langle{A_\tau^2 v_{n+1}},{A_\tau v_{n+1}}\right\rangle \geq 0$$. On the other hand, we have   ⟨Aτvn+τAτwn,vn+τwn⟩=⟨Aτvn,vn⟩+2τ⟨Aτvn,wn⟩+τ2⟨Aτwn,wn⟩. But using an integration by parts in the unbounded part of $$A_\tau$$, we have for all $$\tau \leq \tau_0$$ sufficiently small, and for any $$v,w$$ in $$H^1$$,   |⟨Aτv,w⟩|≤C‖v‖H1‖w‖H1. Hence we have, using the estimates (4.7) on $$w_n$$ (and Definition 2.1),   |⟨Aτvn,wn⟩|≤C⟨Aτvn,vn⟩3/2+C⟨Aτvn,vn⟩1/2⟨Aτvn+1,vn+1⟩ and   |⟨Aτwn,wn⟩|≤C⟨Aτvn,vn⟩2+C⟨Aτvn+1,vn+1⟩2. Gathering together the previous inequalities, we obtain that the sequence $$V_n = \left\langle{A_\tau v_n},{v_n}\right\rangle$$ satisfies the induction formula   (1+2cτ)Vn+1≤Vn+CτVn1/2Vn+1+CτVn3/2+Cτ2Vn2+Cτ2Vn+12 (4.8) and $$V_0 = M^2$$ (after a slight modification of the constant $$M$$). We will prove by induction that we can choose $$a$$ such that for all $$n \in {\mathbb{N}}$$,   Vn=⟨Aτvn,vn⟩≤M2e−anτ. Let us assume that this inequality holds for a given $$n$$. Using (4.4), we have that $$V_{n+1} \leq (C M)^2$$ for some constant $$C$$ independent of $$M$$ and $$\tau \leq \tau_0$$. Hence equation (4.8) implies that   (1+2cτ)Vn+1 ≤Vn+CτMVn+1+CτMVn+Cτ2M2Vn+Cτ2M2Vn+1 ≤ (1+CτM)Vn+CτMVn+1 after possible modifications of the constant $$C$$, and provided $$M$$ and $$\tau_0$$ are small enough. The previous inequality implies   Vn+1≤1+CMτ1+cτ−CMτVn+1. If $$M \tau_0$$ is small enough, we have for all $$\tau \leq \tau_0$$,   1+CMτ1+2cτ−CMτ≤e−aτ for some constant $$a$$ independent of $$\tau$$. This implies that $$V_{n+1} \leq {\it \,{e}}^{-a} V_n \leq M^2 e^{- a (n+1)}$$ and this concludes the induction argument. □ Remark 4.4 If we choose another discretization—implicit–explicit or fully implicit; see (1.7)—then as explained in Remark 4.2 the term in $$\eta^3$$ does not disappear in the equation, and after some manipulation we end up with a recursion formula of the form   vn+1=vn+τ2g−τAτvn+1+τρn with $$\rho_n = \mathcal{O}( \left\|{v_n,v_{n+1}}\right\|_{H^1}^2 ) $$ and $$g \in W$$ a nonzero function. It is relatively easy to prove that in this case, $$v_n$$ converges toward a function $$v_\tau^\infty = \tau A_{\tau}^{-1} g + \mathcal{O}(\tau^2)$$. Indeed, the previous equation is a discretization of the modified problem   ∂tv=τg−Aτv+B(v) whose solution exponentially converges toward the solution of the problem   τg−Aτv+B(v)=0, which exists by the standard implicit function theorem in $$H^1$$. This shows that in this case, $$\psi_n$$ converges exponentially in $$H^1$$ toward a modified soliton $$\eta_{\tau} = \eta + \mathcal{O}(\tau)$$. As the result is weaker than for the linearly implicit scheme, we do not give more mathematical details. 5. Fully discrete case The goal of this last section is to show that the previous results carry over to fully discrete systems. This fact relies essentially on the results in Bambusi & Penati (2010) and Bambusi et al. (2013). After discretization of the previous system by the finite difference method, there exists a discrete ground state $$\eta_h$$ minimizing a discrete convex functional $${\mathcal{H}}_h$$, which is an approximation of the exact functional $${\mathcal{H}}$$ defined above. Here, $$h$$ denotes the space discretization parameter and $$\eta_h$$ is close to $$\eta$$ with an error of order $$h$$. Moreover, the same holds true after Dirichlet cut-off upon adding an exponentially decreasing error term, and it can be proven that the discrete functionals satisfy the same estimates as the continuous ones, uniformly with respect to the discretization parameters. Having fixed a positive parameter $$h$$ we discretize in space by substituting the sequence $$\psi^{{\ell}}\simeq\psi(h {\ell})$$, $${\ell}\in{\mathbb{Z}}$$ for the function $$\psi(x)$$ and the second-order finite difference operator $${\it{\Delta}}_h$$ defined by   (Δhψ)ℓ:=ψℓ+1+ψℓ−1−2ψℓh2 (5.1) for the Laplace operator $$\partial_{xx}$$. We also impose Dirichlet boundary conditions for $$|j| \geq K+1$$, and the parabolic equation (1.4) then becomes   {ddtψℓ=12h2(ψℓ+1+ψℓ−1−2ψℓ)+|ψℓ|2ψℓ,ℓ∈Z,|ℓ|≤K,ψℓ=0,|ℓ|≥K+1,  (5.2) where $$t \mapsto \psi(t) = (\psi^\ell(t))_{\ell \in {\mathbb{Z}}}$$ is an application from $${\mathbb{R}}$$ to $${\mathbb{R}}^{\mathbb{Z}}$$. With this equation is associated a Hamiltonian function and a discrete $$L^2$$ norm given by   Hh(ψ)=h∑j∈Z[|ψj−ψj−1h|2−|ψj|42]andNh(ψ)=h∑j∈Z|ψj|2. (5.3) The discrete space of functions is   Vh,K={ψj∈CZ|ψj=ψ−j,ψj=0forj>K}, equipped with the discrete $$H^1$$ norm   ‖ψ‖h2=2h∑j∈Z|ψj+1−ψj|2h2+h∑j∈Z|ψj|2, (5.4) and we also define the scalar product   ⟨ψ,φ⟩h:=h∑j∈Zψjφj. Following Bambusi & Penati (2010) and Bambusi et al. (2013), we identify $$V_h$$ with a finite element subspace of $$H^1({\mathbb{R}})$$. More precisely, defining the function $$s:{\mathbb{R}} \to {\mathbb{R}}$$ by   s(x)={0if|x|>1,x+1if−1≤x≤0,−x+1if0≤x≤1,  (5.5) the identification is done through the map $$i_h: V_h \to H^1({\mathbb{R}})$$ defined by   {ψj}j∈Z↦(ihψ)(x):=∑j∈Zψjs(xh−j). (5.6) The main results in Bambusi & Penati (2010) and Bambusi et al. (2013) can expressed as follows. Theorem 5.1 For $$h$$ sufficiently small and $$K$$ sufficiently large, there exists a unique real minimizer $$\eta_{h,K} = (\eta_{h,k}^j)_{j \in {\mathbb{Z}}} \in V_{h,K}$$ of the functional $$H_h(\psi)$$ over the set   {ψ=(ψj)j∈Z∈Vh,K,Nh(ψ)=1} equipped with the norm (5.4). Moreover, $$\eta_{h,K}$$ satisfies the equation   12(Δhηh,K)ℓ+|ηh,Kℓ|2ηh,Kℓ =λhηh,Kℓ,ℓ=−K,…,K,ηh,Kℓ =0,|ℓ|≥K+1 (5.7) for some $$\lambda_h > 0$$. Moreover, we can define a local coordinate system $$\psi = \chi_h(r,u) = (1 + r) \eta_{h,K} + u$$ in a vicinity of $$\eta_{h,K}$$ in $$V_{h,K}$$ with $$u \in W_h := \{\hat{\mathrm{E}}\, u \in V_{h,K},\quad \left\langle{u},{\eta_{h,K}}\right\rangle_h = 0\}$$, and such that if $$r_h(u)$$ is defined by the implicit relation $$N_h( \chi_h(r_h(u),u))) = 1$$ then the functional $${\mathcal{H}}_{h}(u) = H_h(\chi_h(r_h(u),u)))$$ has a unique nondegenerate minimum in $$u = 0$$ in $$V_{h,K}$$ and satisfies a convexity estimate of the form (2.7) uniformly in $$h$$ and $$K$$. Finally, the discrete ground state $$\eta_{h,K}$$ is close to the continuous one in the sense that   ‖ihηh,K−η‖H1≤C(h+1h2e−C1hK), (5.8) where $$C$$ and $$C_1$$ do not depend on $$h$$ and where $$i_h$$ is the embedding (5.6). The fully discrete algorithm corresponding to (1.5) and (1.6) then consists of the following two steps. From $$(\psi^\ell_n)_{j = -K}^K$$, such that $$N_h(\psi_n) = 1$$, (i) compute $$\psi^*_n = (\psi^{*,\ell}_n)_{j \in {\mathbb{Z}}}$$ defined by the relation   {ψn∗,ℓ=ψnℓ+τ(12(Δhψn∗)ℓ+|ψnℓ|2ψn∗,ℓ),ℓ∈Z,|ℓ|≤K,ψ∗,ℓ=0,|ℓ|≥K+1;  (5.9) (ii) normalize:   ψn+1ℓ=ψn∗,ℓNh(ψn∗). (5.10) Equation (5.7) guarantees that this fully discrete algorithm preserves exactly the discrete ground state $$\eta_{h,K}$$. We can now copy the proof of the continuous case and adapt it directly to the fully discrete case: we can prove that   ‖ψn−ηh,K‖h≤Ce−cnτ,n≥0, for some constants $$C$$ and $$c$$ independent of $$h$$, $$K$$ and $$\tau$$. Using the previous estimate, we obtain the following result. Theorem 5.2 There exist constants $$M$$, $$C$$, $$C_1$$ and $$c$$ such that for $$h$$, $$\tau$$ sufficiently small and $$K$$ sufficiently large, if $$(\psi_0) = (\psi^\ell_0)_{\ell = -K}^{K}$$ satisfies   ‖ihψ0−η‖H1≤M then the solution $$\psi_n$$ of the fully discrete algorithm (5.9) and (5.10) satisfies   ‖ihψn−η‖H1≤C(e−cnτ+h+1h2e−C1hK). Remark 5.3 Following Remark 4.4, we can prove a similar result for the implicit–explicit and fully implicit algorithms (1.7), but the fully discrete schemes will converge toward a discrete ground state $$\eta_{\tau,h,K}$$ satisfying   ‖ihητ,h,K−η‖H1≤C(τ+h+1h2e−C1hK), and the previous result has to be modified accordingly in these cases. Funding ERC Starting Grant GEOPARDI (No. 279389). References Adhikari S. K. ( 2000a) Numerical solution of the two-dimensional Gross-Pitaevskii equation for trapped interacting atoms. Phys. Lett. A , 265, 91– 96. Google Scholar CrossRef Search ADS   Adhikari S. K. ( 2000b) Numerical solution of spherically symmetric Gross-Pitaevskii equation in two space dimensions. Phys. Lett. E , 62, 2937– 2944. Bambusi D., Faou E. & Grébert B. ( 2013) Existence and stability of ground states for fully discrete approximations of the nonlinear Schrödinger equation. Numer. Math. , 123, 461– 492. Google Scholar CrossRef Search ADS   Bambusi D. & Penati T. ( 2010) Continuous approximation of breathers in one and two dimensional DNLS lattices. Nonlinearity , 23, 143– 157. Google Scholar CrossRef Search ADS   Bao W. & Du Q. ( 2004) Computing the ground state solution of Bose-Einstein condensates by a normalized gradient flow. SIAM J. Sci. Comput. , 25, 1674– 1697. Google Scholar CrossRef Search ADS   Bao W. & Tang W. ( 2003) Ground state solution of trapped interacting Bose-Einstein condensate by directly minimizing the energy functional. J. Comput. Phys. , 187, 230– 254. Google Scholar CrossRef Search ADS   Chiofalo M. L., Succi S. & Tosi M. P. ( 2000) Ground state of trapped interacting Bose-Einstein condensates by an explicit imaginary time algorithm. Phys. Rev. E , 62, 7438– 7444. Google Scholar CrossRef Search ADS   Edwards M. & Burnett K. ( 1995) Numerical solution of the nonlinear Schrödinger equation for small samples of trapped neutral atoms. Phys. Rev. , 1, 1382– 1386. Google Scholar CrossRef Search ADS   Froehlich J., Gustafson S., Jonsson L. & Sigal I. M. ( 2004) Solitary wave dynamics in an external potential. Comm. Math. Phys. , 250, 613– 642. Google Scholar CrossRef Search ADS   Grillakis M., Shatah H. & Strauss W. ( 1987) Stability theory of solitary waves in the presence of symmetry. I. J. Funct. Anal. , 74, 160– 197. Google Scholar CrossRef Search ADS   Grillakis M., Shatah H. & Strauss W. ( 1990) Stability theory of solitary waves in the presence of symmetry. II. J. Funct. Anal. , 94, 308– 348. Google Scholar CrossRef Search ADS   Hairer E. & Lubich C. ( 2014) Energy-diminishing integration of gradient systems. IMA J. Numer. Anal. , 34, 452– 461. Google Scholar CrossRef Search ADS   Weinstein M. I. ( 1985) Modulational stability of ground states of nonlinear Schrödinger equations. SIAM J. Math. Anal. , 16, 472– 491. Google Scholar CrossRef Search ADS   © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.

Journal

IMA Journal of Numerical Analysis – Oxford University Press

Published: Jan 1, 2018

There are no references for this article.