Numerical approximation of the Schr\"odinger equation with concentrated potential

TL;DR

This paper introduces a new family of efficient algorithms for numerically approximating the Schrödinger equation with concentrated potentials, significantly reducing computational costs and enabling longer, more reliable simulations of highly oscillatory solutions.

Contribution

The paper develops tailored convolution quadrature algorithms for Schrödinger equations with concentrated potentials, improving efficiency and accuracy over existing methods.

Findings

Reduced computational cost and memory usage compared to existing approaches.

Ability to perform longer simulations of highly oscillatory solutions.

Successful numerical experiments demonstrating the algorithm's effectiveness.

Abstract

We present a family of algorithms for the numerical approximation of the Schr\"odinger equation with potential concentrated at a finite set of points. Our methods belong to the so-called fast and oblivious convolution quadrature algorithms. These algorithms are special implementations of Lubich's Convolution Quadrature which allow, for certain applications in particular parabolic problems, to significantly reduce the computational cost and memory requirements. Recently it has been noticed that their use can be extended to some hyperbolic problems. Here we propose a new family of such efficient algorithms tailored to the features of the Green's function for Schr\"odinger equations. In this way, we are able to keep the computational cost and the storage requirements significantly below existing approaches. These features allow us to perform reliable numerical simulations for longer times even in cases where the solution becomes highly oscillatory or seems to develop finite time blow-up. We illustrate our new algorithm with several numerical experiments.

Tables

Table 1. Table 1: Error e h subscript 𝑒 ℎ e_{h} for decreasing values of h ℎ h for the linear beating motion problem and the estimated order of convergence.

$h$	2	1	1/2	1/4	1/8
$e^h$	$7.22\times {10}^{-5}$	$3.11\times {10}^{-5}$	$1.12\times {10}^{-5}$	$5.69\times {10}^{-6}$	$2.29\times {10}^{-6}$
EOC		1.2	1.5	0.98	1.3

Equations

$$\frac{1}{\mathrm{i}}{\psi }_t(t,x)=\left(\Delta -\sum _{j=1}^M{V}_j\delta (x-x_j)\right)\psi (t,x),x\in {\mathbb{R}}^D,$$

$$\frac{1}{\mathrm{i}}{\psi }_t(t,x)-\Delta \psi (t,x)=0,x\in {\mathbb{R}}^D,$$

$$k(t,x)=\mathrm{i}{\left(\frac{1}{4\pi \mathrm{i}t}\right)}^{D/2}{e}^{i\Vert x{\Vert }^2/4t},x\in {\mathbb{R}}^D,t>0.$$

$$K(z,x)={\int }_0^{\infty }{e}^{-zt}k(t,x)dt,\mathrm{Re}z>0,$$

$$K(z,x)=\left\{\begin{array}[]{ll}\displaystyle\frac{1}{2\sqrt{z/i}}\exp\left(-|x|\sqrt{z/i}\right)&\text{for }D=1,\\[10.00002pt] \displaystyle\frac{1}{2\pi}K_{0}\left(\|x\|\sqrt{z/i}\right)&\text{for }D=2,\\[10.00002pt] \displaystyle\frac{1}{4\pi\|x\|}\exp\left(-\|x\|\sqrt{z/i}\right)&\text{for }D=3.\end{array}\right.$$

$$\sqrt{z/i}=e^{-i\pi /4}\sqrt{z},$$

$$\psi (t,x)={\int }_{-\infty }^{t}k(t-\tau ,x-x_0)f(\tau ,x_0)d\tau ,$$

$$\frac{1}{\mathrm{i}}{\psi }_t(t,x)-\Delta \psi (t,x)=f(t,x_0).$$

$$u(t)={\int }_0^{t}k(t-\tau )g(\tau )d\tau ,$$

$$|K(z)|\le C|z{|}^{\mu },\text{mbox}for\mathrm{Re}z>0.$$

$$r(z)=1+z{\mathbf{b}}^T(\mathbf{I}-z\mathbf{A}{)}^{-1}\text{mathbbm}1,$$

$$\text{mathbbm}1=(1,1,\dots ,1{)}^T.$$

$${\mathbf{b}}^T{\mathbf{A}}^{-1}=(0,0,\dots ,1).$$

$$\underset{|z|\to \infty }{\lim }r(z)=1-b^T{A}^{-1}\text{mathbbm}1=0$$

$$r(z)=O(z^{-1}),|z|\to \infty .$$

$$\sum _{n=0}^{\infty }{\mathbf{W}}_n{\zeta }^n=K\left(\frac{\mathbf{\Delta }(\zeta )}{h}\right),$$

$$\boldsymbol{\Delta}(\zeta)=\Bigl{(}\mathbf{A}+{\zeta\over 1-\zeta}\mathbbm{1}{\mathbf{b}}^{T}\Bigr{)}^{-1}.$$

$$u_{n+1}=\sum _{j=0}^n\sum _{i=1}^s{\omega }_{n-j}^{i}g(t_j+c_{i}h)=\sum _{j=0}^n{\mathbf{\omega }}_{n-j}{\mathbf{g}}_j,$$

$$(\Delta (\zeta )-zI{)}^{-1}=\sum _{n=0}^{\infty }{\mathbf{E}}_n(z){\zeta }^n,$$

$${\mathbf{\omega }}_n=\frac{h}{2\pi \mathrm{i}}{\int }_{\Gamma }K(z){\mathbf{e}}_n(hz)dz.$$

$${\mathbf{e}}_n(z)=r(z{)}^n\mathbf{q}(z),$$

$${\mathbf{W}}_n=\frac{h}{2\pi \mathrm{i}}{\int }_{\Gamma }K(z){\mathbf{E}}_n(hz)dz.$$

$${\mathbf{E}}_n(z)=r(z{)}^{n-1}(I-z\mathbf{A}{)}^{-1}\text{mathbbm}1{\mathbf{b}}^T(I-z\mathbf{A}{)}^{-1}.$$

$$y^{\prime }(t)=zy(t)+g(t),y(0)=0,$$

$$y_{n+1}(z)=h\sum _{j=0}^n{\mathbf{e}}_{n-j}(hz){\mathbf{g}}_j$$

$$u_{n+1}=\frac{1}{2\pi \mathrm{i}}{\int }_{\Gamma }K(z)y_{n+1}(z)dz.$$

$$\gamma (\xi )=\underset{-\xi \le \mathrm{Re}z\le 0}{\inf }\frac{\log |r(z)|}{\mathrm{Re}z}.$$

$$|r(z)|\le e^{\gamma (\xi )\mathrm{Re}z},$$

$$|r(z)|\le e^{\nu \mathrm{Re}z},\text{for }0\le \mathrm{Re}z\le b,$$

$$\Vert \mathbf{q}(z)\Vert \le C_{\mathbf{q}},\text{for }\mathrm{Re}z\le b,$$

Full text

Keywords: fast and oblivious algorithms, convolution quadrature, Schrödinger equation, boundary integral equations, contour integral methods.

**AMS Classification: ** 65R20, 65L06, 65M15, 65M38

1 Introduction

We consider the efficient numerical approximation of Schrödinger equations with the potential concentrated at a finite set of points in dimension $D=1,2,3$. These problems can be formally described by the equation

[TABLE]

where $\Delta$ denotes the Laplacian and the coupling factors $V_{j}$ may depend on $t$, $V_{j}=V_{j}(t)$, and/or on the value of $\psi$ at $(t,x_{j})$, $V_{j}=V_{j}(\psi(t,x_{j}))$, $j=1,\dots,M$. These models have been used to describe different phenomena in solid state physics, optics and acoustics, and have been rigorously analyzed by several authors in the mathematical physics community, starting from the one-dimensional case [3], followed by the three-dimensional case [1, 2], and more recently the two-dimensional case [11, 12]. The reformulation of these models as $M$-dimensional systems of Volterra integral equations has proven to be very useful for the analysis of the existence, uniqueness, and regularity of solutions. In this paper we address the efficient numerical approximation of such integral representations. Even with the reduction to a finite dimensional system of integral equations, the numerical approximation of (1.1) can be quite delicate [13], particularly in the nonlinear case. We notice that the long-time behaviour of the solution to the nonlinear Schrödinger equation (NLS) with concentrated potential is not always well-understood and it might become highly oscillatory and even develop blow-up in finite time.

Let us recall that for the Schrödinger equation in the whole space

[TABLE]

with $D=1,2,3$, the Green’s function is given by [15, Equation (2.24)]

[TABLE]

The corresponding transfer operator, i.e., the Laplace transform of $k$

[TABLE]

is given by

[TABLE]

In the above $K_{0}$ is the modified Bessel function of the second kind [14] and

[TABLE]

where $\sqrt{z}$ is the branch with the positive real part and with the branch cut along the negative real axis. Hence the transfer operator $K(x,z)$ as a function of $z$ is analytically extended to the cut complex plane $\mathbb{C}\setminus(-\infty,0]$.

In applications, see [19, 13] and Section 6, it is important to be able to compute accurately and efficiently the convolution in time with the Green’s function. Namely,

[TABLE]

solves the Schrödinger equation with a source at $x_{0}$:

[TABLE]

Due to the kernel being non-local and highly oscillatory, accurate and efficient discretization of (1.6) is not easy. In [13] a carefully constructed, accurate numerical method is presented and in [21] numerical experiments illustrated the good computational properties of convolution quadrature (CQ) for (1.6). Due to the non-locality of the kernel, computing $\psi(x,t)$ at $N$ time steps $t_{j}$ using the method in [13] has an $O(N^{2})$ computational complexity, whereas standard FFT based methods for CQ, [20, 10] can reduce this to $O(N\log N)$ for linear problems and to $O(N\log^{2}N)$ for nonlinear problems [16, 10]. However, both methods require to store $N$ solution vectors in memory. In this work we describe an algorithm that is both very easy to implement and can significantly reduce the amount of memory used. More precisely the memory requirements will be reduced from $O(N)$ to $O(n_{0}+\log N)$, with $n_{0}\ll N$. Our method belongs to the family of oblivious algorithms [24, 7], the name indicating that the memory requirements can be significantly reduced. Due to the way our algorithm is built we also expect it to be extendable to a variable step implementation, something which by construction is difficult for FFT based methods [16, 10].

Oblivious quadratures were first developed for parabolic problems [23, 24], where the transfer operator is sectorial, i.e., it admits a holomorphic extension to the complement of an acute sector in the left half of the complex plane and it grows at most algebraically as $|z|\to\infty$. The extension of such ideas to hyperbolic problems, where the transfer operator typically exhibits exponential growth as $\operatorname*{Re}z\to-\infty$, is much more recent and has been first developed for the two-dimensional and the damped three dimensional wave equation [7]. The application to the Schrödinger equation has never been addressed to our knowledge. Furthermore in this paper we use a new approach which allows a substantial simplification of the implementation with respect to the algorithms in [23, 24, 7] and in our experience even a marginal improvement in the compressibility and the computational times. In particular, we take the real inverse Laplace transform approach introduced in [6] for the fractional integral, which fits the more favourable sectorial framework. We thus generalize the ideas in [6] to the non-sectorial situation of the Schrödinger Green’s kernel and propose a mosaic-free, fast and oblivious algorithm for the approximation of (1.6). Here mosaic refers to the special partition of the integration domain $\tau\leq t$ and organization of the book-keeping which is required to compress the memory by the algorithms in [23, 24, 7]. The algorithm we propose here will only require the computation and storage of certain quantities at the beginning of the integration procedure and the update of a unique set of ordinary differential equations from one step to the next one, for the whole time interval. More details about the implementation are given in Section 5.

The outline of the paper is as follows. In Section 2 we give a brief introduction to Runge-Kutta based convolution quadrature, in Section 3 we describe a new (real) integral representation of the convolution quadrature weights associated to the transfer operator (1.4). An efficient quadrature of this representation of the convolution weights is described in Section 4. The use of this quadrature in an efficient algorithm for computing discrete convolutions is explained in Section 5. The new method is then illustrated by several substantial numerical experiments in Section 6. In particular we describe in detail an application to a nonlinear Schrödinger equation describing the suppression of quantum beating taken from [13]. The codes used for these experiments are published at [5].

2 Runge-Kutta convolution quadrature

In this section we briefly describe convolution quadrature (CQ) as applied to the evaluation of one sided convolutions

[TABLE]

where $k$ is a given kernel with Laplace transform $K(z)=\mathscr{L}k(z)$ and given data $g$. In the applications in this paper, $k$ will also be a function of $x$, but as this dependence plays no role in this section we supress it for now. A basic assumption for the application of CQ to (2.1) is that there exist $C>0$ and $\mu\in\mathbb{R}$ such that

[TABLE]

In this paper we use CQ based on implicit $A$-stable Runge-Kutta methods [17]. We employ standard notation for an $s$-stage Runge-Kutta discretization based on the Butcher tableau described by the matrix $\mathbf{A}=(a_{ij})_{i,j=1}^{s}\in\mathbb{R}^{s\times s}$ and the vectors ${\mathbf{b}}=(b_{1},\ldots,b_{s})^{T}\in\mathbb{R}^{s}$ and ${\mathbf{c}}=(c_{1},\ldots,c_{s})^{T}\in[0,1]^{s}$. The corresponding stability function is given by

[TABLE]

where

[TABLE]

Recall that $A$-stability is equivalent to the condition $|r(z)|\leq 1$ for $\operatorname*{Re}z\leq 0$. In the following we collect all the assumptions on the Runge-Kutta method. These are satisfied by, for example, Radau IIA and Lobatto IIIC families of Runge-Kutta methods.

Assumption 1**.**

1. (a)

The Runge-Kutta method is $A$-stable with (classical) order $p\geq 1$ and stage order $q\leq p$. 2. (b)

The stability function satisfies $|r({\mathrm{i}\mkern 1.0mu}y)|<1$ for all real $y\neq 0$. 3. (c)

The Runge-Kutta coefficient matrix $\mathbf{A}$ is invertible. 4. (d)

The Runge-Kutta method is stiffly accurate, i.e.,

[TABLE]

This implies that

[TABLE]

and $c_{s}=1$.

Since $r(z)$ is a rational function, the above assumptions imply that

[TABLE]

Following the theory in [22], we define the weight matrices ${\mathbf{W}}_{n}$ corresponding to the operator $K$ as the coefficients of the power expansion

[TABLE]

where $h$ is the step size and the matrix-valued function $\boldsymbol{\Delta}(\zeta)$ is the so-called symbol of the Runge–Kutta method:

[TABLE]

Denoting by $\boldsymbol{\omega}_{n}=(\omega_{n}^{1},\dots,\omega_{n}^{s})$ the last row of ${\mathbf{W}}_{n}$, the approximation to the convolution integral (2.1) at time $t_{n+1}=(n+1)h$ is given by

[TABLE]

with the column vector {\mathbf{g}}_{j}=g(t_{j}+\mathbf{c}h)=\bigl{(}g(t_{j}+c_{i}h)\bigr{)}_{i=1}^{s}.

The convergence order of this approximation has been investigated in [22] for parabolic problems, i.e., for sectorial $K$, and in [8] and [9] for hyperbolic problems, i.e., for non-sectorial operators.

With the row vector ${\mathbf{e}}_{n}(z)=(e_{n}^{1}(z),\dots,e_{n}^{s}(z))$ defined as the last row of the $s\times s$ matrix ${\mathbf{E}}_{n}(z)$ given by

[TABLE]

we obtain an integral formula for the weights

[TABLE]

This representation follows from Cauchy’s formula and the definition of the weights in (2.5), with the integration contour $\Gamma$ chosen so that it surrounds the poles of ${\mathbf{e}}_{n}(hz)$. An explicit expressing for ${\mathbf{e}}_{n}$ is given by

[TABLE]

with the row vector $\mathbf{q}(z)={\mathbf{b}}^{T}(I-z\mathbf{A})^{-1}$; cf. [22, Lemma 2.4]. The $A$-stability assumption implies that the poles of $r(z)$ are all in the right-half plane. Further, due to the decay of the rational function $r(z)$, see (2.4), for $n>\mu+1$ with $\mu$ in (2.2), the contour $\Gamma$ can be deformed into the imaginary axis.

For the weight matrices it holds

[TABLE]

By [22, Lemma 2.4], for $n\geq 1$, ${\mathbf{E}}_{n}(z)$ is the rank-1 matrix given by

[TABLE]

The Runge-Kutta approximation of the inhomogeneous linear problem

[TABLE]

at time $t_{n+1}$ is given by

[TABLE]

and thus the approximation of the convolution integral in (2.7) can be rewritten as [22, Proposition 2.4]

[TABLE]

We will require the following technical lemmas proved in [7] and [6] where examples of numerically computed values of $\gamma$ can also be found.

Lemma 2**.**

Let $r(z)$ be the stability function of a Runge-Kutta method satisfying Assumption 1 and let

[TABLE]

Then $\gamma(\xi)\in(0,1]$ for $\xi>0$, it monotonically increases as $\xi\to 0$ and

[TABLE]

for all $z$ in the strip $-\xi\leq\operatorname*{Re}z\leq 0$.

Lemma 3**.**

There exist constants $\nu>1$, $b>0$ and $C_{{\mathbf{q}}}>0$ such that

[TABLE]

and

[TABLE]

where $C_{{\mathbf{q}}}$ depends on the choice of the norm $\|\cdot\|$.

Lemma 4**.**

There exists a constant $C_{A}>0$ such that

[TABLE]

Proof.

The estimate follows from $r(z)={\mathbf{b}}^{T}{\mathbf{A}}^{-1}(I-z{\mathbf{A}})^{-1}\mathbbm{1}$, $\mathbf{q}(z)={\mathbf{b}}^{T}(I-z{\mathbf{A}})^{-1}$, and the fact that the eigenvalues of ${\mathbf{A}}$ have a strictly positive real part. ∎

Remark 5**.**

For the backward Euler method, $C_{A}=1$. For other RK methods the constant can be estimated numerically. We obtained that for the 2-stage Radau IIA method $C_{A}\approx 2.1213$, for the 3-stage Radau IIA $C_{A}\approx 3.479$, and for the 3-stage Lobatto IIIA $C_{A}\approx 3.6224$.

In the next section we discuss how to approximate the integral in (2.9) by an efficient quadrature rule.

3 Integral representation of the convolution quadrature weights

We will follow the same idea as in the derivation of the real inversion formula for the Laplace transform in [18, Section 10.7], but with ${\mathbf{e}}_{n}(hz)$ in place of $e^{zt}$. For the rest of the paper we will write

[TABLE]

In the following results we will need some properties of the modified Bessel function $K_{0}$. First of all, by $K_{0}$ we mean the principal branch analytic in the cut complex plane $\mathbb{C}\setminus(-\infty,0]$ as described in [14, §10.25]. The large argument behaviour is

[TABLE]

see [14, 10.25.3]. Whereas the small argument behaviour is

[TABLE]

see [14, 10.25.3]. From these two results it follows that there exists a constant $C_{K_{0}}>0$ such that

[TABLE]

Lemma 6**.**

For $n\geq 1$,

[TABLE]

Proof.

Note that for $\operatorname*{Re}z\geq 0$, $\operatorname*{Re}\sqrt{z}e^{-{\mathrm{i}\mkern 1.0mu}\pi/4}\geq 0$ and hence $|K(z,x)|\leq\frac{1}{2}|z|^{-1/2}$ for $D=1$, $|K(z,x)|\leq C_{K_{0}}|z|^{-1/2}$ for $D=2$ and $|K(z,x)|\leq\frac{1}{4\pi d}$ for $D=3$. As $\|{\mathbf{e}}_{n}(z)\|\leq C_{A}^{n+1}|z|^{-n-1}$ we can deform the contour as required by the statement of the lemma. ∎

Lemma 7**.**

Given $d>0$ and $\xi>0$

[TABLE]

are decreasing functions of $y\geq 0$.

Proof.

Note that it is sufficient to show that $g_{\pm}(y)$ decreases with $g_{\pm}$ the real part of the exponent:

[TABLE]

Let us first consider $g_{-}(y)$:

[TABLE]

Notice that the cosine above is negative. Then by using the half-angle formula for the cosine we obtain

[TABLE]

As $g^{\prime}_{-}(y)=-\frac{d}{2\sqrt{2}}\sqrt{(\xi^{2}+y^{2})^{1/2}-y}/\sqrt{\xi^{2}+y^{2}}\leq 0$, the function $f_{-}(y)$ decreases as $y$ increases.

For the exponent of $f_{+}(y)$ we obtain

[TABLE]

Clearly this function is decreasing, hence again $f_{+}(y)$ decreases as $y$ increases. ∎

Theorem 8**.**

The weights are given by

[TABLE]

where

[TABLE]

Further, the function $G(\lambda,x):=K(\lambda e^{-i\pi},x)-K(\lambda e^{i\pi},x)$ is given by

[TABLE]

Proof.

We first deform the integration contour from the imaginary axis to the contour described in Figure 1. Letting $R\rightarrow\infty$ and $\delta\rightarrow 0$ and using the estimate for $K(z)$ from the previous lemma for $D=1$ and $D=2$ and the bound on $e_{n}(z)$ as before we obtain the required expression (3.5). As (3.4) shows that the kernel in $D=2$ dimensions is bounded up to a constant by the same expression as in $D=1$ dimension, the same argument works here too.

To derive the simplified expression for $G(\lambda,x)$ for the 3D Schrödinger equation note that

[TABLE]

Similarly in 1D we have

[TABLE]

In 2D we have

[TABLE]

where in the last step we used [14, 10.27.9]. ∎

4 Quadrature for the convolution weights

In this section we develop an efficient quadrature for the approximation of the convolution quadrature weights, written as (3.5). Let $T$ be the final time, $h=T/N$ the time-step and $N$ the total number of time steps in our approximation (see (2.7)). Our goal is to use the same quadrature weights and nodes for the approximation of $\boldsymbol{\omega}_{n}(x)$ for all $n_{0}<n\leq N$, for some $n_{0}\ll N$ to be determined, and all $0\leq x\leq L$, for $L$ the maximal distance between the points $x_{j}$ in (1.1).

In the first place we will bound the size of the contribution along the vertical lines $\Gamma_{\pm}$, showing that in some cases it can be neglected. Thus we will approximate

[TABLE]

For the rest of the paper we restrict ourselves to $D=1$. In view of the very similar expression for the 3D kernel and the bound (3.4) in 2D, we expect that similar estimates will hold also in higher dimensions. The single bigger change would be the treatment of the singularity at the origin.

4.1 Truncation to a finite interval

We start by bounding the contribution to the integral along the vertical semilines, i.e. the error in (4.1).

Proposition 9**.**

We have the bound

[TABLE]

or more explicitly

[TABLE]

Proof.

To prove the result we need to bound the integrals over $\Gamma_{+}\cup\Gamma_{-}$ in Theorem 8. The estimate (4.2) follows directly from the definition of the integrand.

Note that Lemma 7 implies

[TABLE]

and

[TABLE]

We require a bound on ${\mathbf{e}}_{n}(h(-\xi\pm{\mathrm{i}\mkern 1.0mu}y))$ which follows from ${\mathbf{e}}_{n}(z)=r(z)^{n}{\mathbf{q}}(z)$, Lemma 2 and Lemma 4

[TABLE]

Hence

[TABLE]

Similarly

[TABLE]

∎

The following corollary will be used to determine the value of $n_{0}$. More details are given in Section 5.

Corollary 10**.**

For a given $\varepsilon>0$ and $\xi>1$

[TABLE]

if

[TABLE]

In order to efficiently approximate ${\mathbf{I}}_{n}$, we use the splitting

[TABLE]

with

[TABLE]

where $L_{-1}=0$, $L_{0}>0$ is a free parameter, and $L_{j}=(1+B)L_{j-1}$ for $j>1$ and some fixed $B\geq 1$. Every sub-integral ${\mathbf{I}}_{n,j}$ will be approximated by an appropriate Gauss quadrature. The case of ${\mathbf{I}}_{n,0}$ is treated separately due to the integrable singularity of the integrand at 0. To analyse the error due to Gauss quadrature, we use the following classical result.

Theorem 11**.**

Let $f$ be analytic inside the Bernstein ellipse

[TABLE]

with $\varrho>1$ and bounded there by $M$. Then the error of Gauss quadrature with weight $w(x)$ is bounded by

[TABLE]

where $If=\int_{-1}^{1}w(x)f(x)dx$ and $I_{Q}f=\sum_{j=1}^{Q}w_{j}f(x_{j})$ is the corresponding Gauss formula, with weights $w_{j}>0$.

Proof.

To the best of our knowledge, the first proof of this result appeared in [26]. A different proof for $w(x)\equiv 1$, which can easily be extended to the case of a general weight as in [6], can be found in [25, Chapter 19]. ∎

4.2 Gauss-Jacobi quadrature for the initial interval

We fix the first interval $[0,L_{0}]$ and compute

[TABLE]

In the 1D case, $G(\lambda,x)\sqrt{\lambda}={\mathrm{i}\mkern 1.0mu}e^{{\mathrm{i}\mkern 1.0mu}\pi/4}\cosh(de^{{\mathrm{i}\mkern 1.0mu}\pi/4}\sqrt{\lambda})$ is an entire function of $\lambda$. Hence in this interval we will use Gauss-Jacobi quadrature with weight $w(x)=(x+1)^{-1/2}$ on the interval $[-1,1]$. We denote by $\boldsymbol{\tau}_{n,0}(Q)$ the corresponding quadrature error when taking $Q$ quadrature nodes.

To estimate the error of the quadrature, we will need to bound, according to Theorem 11,

[TABLE]

where we have already neglected the modulus one quantity ${\operatorname*{e}}^{{\mathrm{i}\mkern 1.0mu}\pi/4}$ included in the definition of $G(\lambda,d)$. Notice that there is a maximal value of $\varrho$, which we denote $\varrho_{\max}$, determined by the location of the poles of ${\mathbf{e}}_{n}$.

Theorem 12**.**

Let $b$ and $\nu$ be as in Lemma 3,

[TABLE]

and

[TABLE]

Then, if $\varrho_{opt}\in(2+\sqrt{3},\varrho_{\max})$, we can bound the error of the $Q$-node Gauss-Jacobi quadrature by

[TABLE]

Otherwise

[TABLE]

Remark 13**.**

The second estimate, seems quite pessimistic because of the exponentially growing term. However, notice that we can take $b=h$ which implies $\nu\approx 1$ and $\varrho_{\max}>1+4/L_{0}$. Hence both estimates imply that if we choose $L_{0}\propto T^{-1}$, $Q=O(\log\varepsilon)$ quadrature nodes are sufficient to obtain $\|\boldsymbol{\tau}_{n,0}(Q)\|\leq\varepsilon$. Note, that in numerical experiments reported in this paper, it was always the case that $\varrho_{opt}\in(2+\sqrt{3},\varrho_{\max})$.

Theorem 12.

Note that

[TABLE]

As we will want to avoid the poles of ${\mathbf{e}}_{n}$ and use Lemma 3 we need that $\varrho\leq\varrho_{\max}$. This upper bound (4.7) is obtained as a solution of

[TABLE]

From Lemma 3 and the definition of ${\mathbf{e}}_{n}$ we can now bound

[TABLE]

Let $\varrho=e^{\delta}$ for $\delta>0$. Then

[TABLE]

Assuming $\varrho\geq 2+\sqrt{3}$ implies $\cosh\delta\geq 2$ and

[TABLE]

For $\operatorname*{Re}\zeta\geq-1,\zeta\in\mathcal{E}_{\varrho}$ we have that $|r(-h(1+\zeta)L_{0}/2)|\leq 1$ and hence

[TABLE]

where we have used that $(1+x)^{1/2}\leq\sqrt{3}(x-1)$ for $x\geq 2$.

Therefore, from Theorem 11 we deduce that

[TABLE]

So we minimize

[TABLE]

As

[TABLE]

and

[TABLE]

the minimum is reached at

[TABLE]

Using the identities

[TABLE]

we obtain the value of $\varrho_{opt}={\operatorname*{e}}^{\delta_{opt}}$ in the statement and

[TABLE]

Using now that $-1+\sqrt{1+x^{2}}\leq x$, for $x\geq 0$, we have

[TABLE]

Hence, in case $\varrho_{opt}\in(1,\varrho_{\max})$ the following bound holds

[TABLE]

Otherwise we choose $\varrho=\varrho_{\max}={\operatorname*{e}}^{\delta_{\max}}$ and obtain

[TABLE]

and the stated bound for the error. ∎

4.3 Gauss quadrature away from the singularity

In this section we analyze the error in the Gauss-Legendre ($w(x)=1$) quadrature of the integrals ${\mathbf{I}}_{n,j}$ in (4.4), with $j\geq 1$, which can be written as

[TABLE]

with $\Delta L_{j}:=L_{j}-L_{j-1}=BL_{j-1}$.

Theorem 14**.**

Let $\boldsymbol{\tau}_{n,j}(Q)$ be the error in the approximation of ${\mathbf{I}}_{n,j}$ by Gauss quadrature with weight $w(x)=1$ and $Q$ quadrature nodes. Then

[TABLE]

with

[TABLE]

and

[TABLE]

and

[TABLE]

Proof.

According to Theorem 11 we now need to bound the function

[TABLE]

In order to avoid the singularity of the square root we require

[TABLE]

which is satisfied for $1<\varrho<\varrho_{\max}$ and

[TABLE]

Note that (using $|z|\leq\sqrt{1+(\operatorname*{Im}z/\operatorname*{Re}z)^{2}}|\operatorname*{Re}z|$)

[TABLE]

where $\eta_{\theta}=1+B((\varrho+\varrho^{-1})\cos\theta+2)/4$.

With this notation we can bound, for every $\zeta\in\mathcal{E}_{\varrho}$, $\varrho\in(1,\varrho_{\max})$

[TABLE]

The result then follows from Theorem 11. ∎

Remark 15**.**

The above result is somewhat unsatisfactory as it still contains a min-max problem. The following corollary simplifies the estimate but is too pessimistic in practice. Hence, we make use of the corollary only for the discussion about the complexity of the algorithm and in practice numerically solve the above min-max problem in order to obtain optimal parameters.

Corollary 16**.**

With notation as in Theorem 14

[TABLE]

where

[TABLE]

Proof.

The proof follows from the fact that $\gamma(\xi)$ decreases for increasing argument and that $\eta_{\theta}$ from Theorem 14 decreases from $\theta=0$ to $\theta=\pi$. ∎

Remark 17**.**

To understand the required number of quadrature points let us set $B=3$, as in the numerical experiments, giving $\varrho_{\max}=3$. Choosing $\varrho=2$ we get $\eta_{+}=35/8$ and $\eta_{-}=5/8$ and

[TABLE]

where $C=\frac{12}{\pi}\eta_{-}^{-1/2}$, $c_{+}=2^{1/4}\eta_{+}^{1/2}$ and $c_{-}=\gamma(hL_{j-1}\eta_{+})\eta_{-}$. Hence to obtain accuracy $\varepsilon$ we need to choose $Q=O(\log\varepsilon^{-1}+\frac{d^{2}}{t_{n}})$ quadrature nodes. Furthermore, this indicates that the integrand can get exponentially large as $t_{n}\rightarrow 0$ which can create difficulties in finite precision arithmetic.

5 Parameter choice and the fast method

In this section we specify the choice of the most relevant parameters in our method and explain our algorithm.

In the first place we choose the truncation parameter $\xi$. The estimate in Corollary 10 suggests choosing $\xi\propto h^{-1}$. Hence, let $\xi=A_{0}/h$ for some constant $A_{0}>0$. Then if $n_{0}$ is chosen such that

[TABLE]

the error due to truncation is less than $\varepsilon_{0}$ for $n\geq n_{0}$. Note that this means that at least $O(h^{-1/2})$ weights will need to be computed directly. Here we use the estimate (4.3). Some minor gains can be made by computing the truncation error numerically using (4.2).

Next, according to Theorem 12 we need to choose $L_{0}$ proportional to $T^{-1}$, for $T$ the final time in our approximation as stated at the beginning of Section 4. Hence we set $L_{0}=A_{1}/T$ for some constant $A_{1}>0$. There will then be $J$ intervals away from the singularity where Gauss quadrature is used with $L_{0}(1+B)^{J}=\xi$ implying $J=\frac{1}{(1+B)}\log\left(\frac{A_{0}T}{h}\right)=O(\log(T/h))$.

Once all these parameters are set we can choose the number of quadrature nodes and weights in each interval so that each quadrature error $\boldsymbol{\tau}_{n,j}(Q)\leq\varepsilon$ for $j=0,\dots,J$ and $n\geq n_{0}$. Then with the choice $\varepsilon_{0}=\mathrm{tol}/2$ and $\varepsilon=\frac{\mathrm{tol}}{2(J+1)}$ we have that the total error is bounded by the tolerance $\mathrm{tol}>0$. As explained in Remark 13, $Q=O(\log\varepsilon)$ guadrature nodes of the Gauss-Jacobi quadrature are sufficient to obtain error of size $\varepsilon$ in this first interval. According to Remark 17, we need to use $O\left(\log\varepsilon^{-1}+\frac{d^{2}}{n_{0}h}\right)$ quadrature nodes in each of the remaining intervals.

Let us now describe the fast method for computing the discrete convolution ${\mathbf{u}}_{n}=\sum_{j=0}^{n}\boldsymbol{\omega}_{j}(x){\mathbf{f}}_{n-j}$ for a general sequence of $s\times 1$ vectors ${\mathbf{f}}_{j}$, $j=0,\dots,N$. Using this formula requires the full sequences $\boldsymbol{\omega}_{j}$ and ${\mathbf{f}}_{j}$, for $j=0,\dots,N$ to be kept in memory requiring $O(N)$ memory. This can instead be evaluated most efficiently by applying our quadrature approximation of the convolution weights in the same way as in [6, Section 7.2]. We thus split the sum into two terms, the local term, with summation index $j=1,\dots,n_{0}$, and the remaining history term:

[TABLE]

The local term is evaluated directly, by precomputing and keeping in memory the first $n_{0}+1$ convolution weights. The history term is instead evaluated by means of a fast summation algorithm which is based on the quadrature developed in Section 4. After replacing the CQ weights by the result of applying our quadrature we are led to the formula

[TABLE]

with

[TABLE]

satisfying the recursion

[TABLE]

Let us investigate the complexity and the memory requirements of the above described algorithm. The vectors ${\mathbf{f}}_{j}$, $j=n-n_{0}-1,n-n_{0},\dots,n$ and ${\mathbf{Q}}_{n-1,k}$, $k=1,\dots,N_{Q}$ need to be kept in memory at each time step $t_{n}$, i.e., $n_{0}+N_{Q}$ vectors need to be stored at any time. Recalling now that $N_{Q}=O\left(\left(\log\varepsilon^{-1}+\frac{d^{2}}{n_{0}h}\right)\log(T/h)\right)$ the total memory requirements are given by $O\left(n_{0}+\left(\log\varepsilon^{-1}+\frac{d^{2}}{n_{0}h}\right)\log(T/h)\right)$. If we choose $n_{0}=O(h^{-1/2})$, then the memory requirements are reduced to $O((h^{-1/2}+\log\varepsilon^{-1})\log h^{-1})$. Turning to computational complexity, if the local term is computed using FFT methods as described in [10, 16], the total computational cost is $O(n_{0}\log^{2}n_{0}+NN_{Q})$. This reduces to $O(h^{-3/2}\log h^{-1})$ if we choose $n_{0}=O(h^{-1/2})$.

We need to say a few more words about the choice of $n_{0}$. Any choice of $n_{0}$ as a function of $h$ that implies $t_{n_{0}}\rightarrow 0$ and $h\rightarrow 0$ would allow for exponentially increasing integrand, see Remark 17. In finite precision arithmetic this could lead to destructive cancellation and complete loss of accuracy. Choosing $n_{0}=O(h^{-1})$ would entirely remove this difficulty however would require us to use many more direct steps and more memory. Nevertheless, even so our algorithm would bring many advantages if $n_{0}\ll N$; compare this with [7]. Hence, $n_{0}$ should be chosen between $\text{\em const}\cdot h^{-1/2}$ and $\text{\em const}\cdot h^{-1}$ depending on the parameters of the problem investigated. In the two applications that we investigate in the numerical experiments, we choose the smallest possible $n_{0}$ that ensures that the truncation error, see Proposition 9, is bounded by the tolerance, i.e., $n_{0}=O(h^{-1/2})$.

6 Numerical experiments

Next, we present results of various experiments to illustrate the new method. The codes used to perform these experiments can be found at [5].

6.1 Approximation of CQ weights

In the first place we test our new quadrature and compare the CQ weights we obtain with those given by the standard method based on FFT.

We start by setting the parameters as in Section 5 with $A_{0}=1$ and $A_{1}=2$. We fix $B=3$ and given a tolerance $\mathrm{tol}>0$ compute the number of nodes in each interval so that the error is bounded as described in Section 5. The results for $h=10^{-2},10^{-3}$, $T=100$, $d=1$ are given in Figure 2. As reported in the caption of Figure 2 and the legends, for a target accuracy of 3 digits in the computation of the CQ weights we need $n_{0}=8$ and $N_{Q}=28$ for $h=10^{-2}$, and $n_{0}=20$, $N_{Q}=48$ for $h=10^{-3}$. Thus the storage in our algorithm will be reduced from $N=10^{4}$ to $8+28$ quantities if $h=10^{-2}$, and from $N=10^{5}$ to $20+48$ if $h=10^{-3}$. For a more stringent target accuracy of six digits in the computation of the CQ weights, we need instead 15+48 quantities in storage (rather than $10^{4}$) for $h=10^{-2}$, and 27+74 for $h=10^{-3}$ (rather than $10^{5}$).

We can see that the chosen tolerance is not exceeded and the error for $n$ close to $n_{0}$ is close to the tolerance. However, the error for larger $n$ is much better than the tolerance. The truncation error and the quadrature error on the intervals away from the singularity all get quickly better with the larger $n$ so this is not so surprising. However, as the error bound for the Gauss-Jacobi quadrature in the first interval does not improve with increasing $n$, this suggests that possibly our estimate in the first interval is not optimal. This is not a great issue, as this single quadrature does not contribute a great deal to the overall costs: around 5% of quadrature points are in this interval in the above calculations.

6.2 A linear Schrödinger equation with concentrated potential

We start by considering the same example as in [21] and compute the solution $\psi(x,t)$, for $x\in\mathbb{R}$, $t>0$ to:

[TABLE]

for some given time-dependent amplitudes $V_{j}(t)$.

Assuming that

[TABLE]

the values $q_{j}(t):=\psi(x_{j},t)$ satisfy the system of Volterra integral equations

[TABLE]

for $j=1,\dots,M$ and the solution to (6.1) can be written

[TABLE]

with $k$ as in (1.3) with $D=1$. For more detail on derivation of this system see [21].

For our experiments we take $M=2$,

[TABLE]

The application of Runge–Kutta based CQ to (6.2) yields for $j=1,\dots,M$ approximations ${\mathbf{q}}_{j,n}\approx q_{j}({\mathbf{t}}_{n})$, with ${\mathbf{t}}_{n}=(t_{n}+c_{i}h)_{i=1}^{s}$, $n=0,1,\dots,N$, defined by

[TABLE]

where ${\mathbf{V}}_{k,\ell}:={\rm diag}(V_{k}(t_{\ell}+c_{1}h),\dots,V_{k}(t_{\ell}+c_{s}h))$.

We then need to solve the discrete linear system

[TABLE]

This is

[TABLE]

To compute the memory term on the right-hand side we use the fast method described in Section 5. The optimal choice of the various Gauss quadrature weights and nodes used in this algorithm depend on $d$ which is above $d=0$ or $d=2$. However from our analysis it follows that the error when using the optimal choice for a certain $d=d_{\max}$ will be bounded by the tolerance for every $d\in(0,d_{\max})$. We thus compute the optimal quadrature for $d=2$ and use it to represent both ${\mathbf{W}}_{j}(0)$ and ${\mathbf{W}}_{j}(2)$, for every $j=n_{0}+1,\dots,N$.

In Figure 3 we show a plot of the solution and the difference in the solutions obtained by our fast method and the standard implementation of CQ. In these experiments we choose $T=40$, $h=0.1$, $\mathrm{tol}=10^{-6}$ and $d=2$. Our optimization routines returned $n_{0}=14$ and the total number of quadrature points $N_{Q}=36$. As we see, the error is slightly larger than the tolerance, which is to be expected as the tolerance is valid for the computation of the weights not the final result of the discrete convolution.

6.3 A non linear Schrödinger equation with concentrated potential

We consider nonlinear Schrödinger equations with concentrated potentials and, in particular, the case studied in [13]. Thus, we compute the solution $\psi(x,t)$, for $x\in\mathbb{R}$, $t>0$ to:

[TABLE]

with $\gamma_{j}<0,\ \sigma\geq 0$. Writing the solution using Duhamel’s principle gives

[TABLE]

where $k$ is the Green’s function (1.3) for $D=1$ and $\phi(t,x)$ is the solution of the homogeneous problem:

[TABLE]

Evaluating (6.8) at $x=x_{k}$, $k=1,\dots,M$, gives the following system of integral equations

[TABLE]

In the following numerical experiments we set $M=2,\ x_{1}=-a,\ x_{2}=a,\ \gamma_{1}=\gamma_{2}=\gamma<0$,

[TABLE]

giving the system of equations

[TABLE]

For optimal performance, convolution quadrature requires that the data can be extended smoothly to negative times by zero. Since $q_{1,2}(0)\neq 0$ we modify the system as follows

[TABLE]

$k=1,2$ with correction terms

[TABLE]

To compute the correction terms we use the formula obtained using symbolic computation software

[TABLE]

where ${\mathrm{erf}\mkern 1.0mu}$ is the error function.

After discretization of (6.13) using CQ as in Section 6.2 we obtain a non-linear system to be solved at each step of the form

[TABLE]

where $\mathbf{H}(t)$ is the history term containing terms known at time-step $n$. We solve the non-linear equation by a fixed-point iteration with the initial guess given by the solution at the previous time-step. The history term $\mathbf{H}({\mathbf{t}}_{n})$ is computed using our fast method.

6.4 The linear case

Let us first consider the linear case, i.e., $\sigma=0$. We look for solutions of the form

[TABLE]

with $\alpha^{2}+\beta^{2}=1$ and

[TABLE]

and

[TABLE]

where constants $N_{\mathrm{f}}$ and $N_{\mathrm{e}}$ are chosen so that $\|\phi_{\mathrm{f}}\|_{L^{2}(\mathbb{R})}=\|\phi_{\mathrm{e}}\|_{L^{2}(\mathbb{R})}=1$. Here, $\phi_{\mathrm{f}}$ and $\phi_{\mathrm{e}}$ correspond to the fundamental respectively excited state, see [13]. The initial data is hence of the form

[TABLE]

Substituting this ansatz into (6.7) with $\sigma=0$, gives the following relations that need to be satisfied by $\lambda_{\mathrm{f}}$ and $\lambda_{\mathrm{e}}$

[TABLE]

and

[TABLE]

We choose the parameters $\alpha=\sqrt{0.01}$, $\beta=\sqrt{0.99}$, $a=3$, $\gamma=-0.5$. Solving numerically the above nonlinear equations gives the eigenvalues $\lambda_{f}=0.085894322668323$ and $\lambda_{e}=0.021229338264198$. For this special initial data, formulas similar to (6.14) are available for expressing the solution to the free Schrödinger equation $\phi(t,x)$ in terms of the error function, see [13].

We apply the fast convolution quadrature method based on the 2-stage Radau IIA Runge-Kutta method. The numerical results with $T=100$ and $h=1$ are given in Figure 4. There we plot $|q_{j}(t)|^{2}$ and the error

[TABLE]

where $q_{j}(t)$ is the numerically obtained solution and $q^{\text{ex}}_{j}(t)$ the exact solution. The results are of high quality with, as expected, larger error near $t=0$ and slight increase in error with increase in time.

The importance of this example is that it shows the so-called beating motion of the system. This is not an easy problem to solve numerically and these excellent results in the linear case give us confidence in the nonlinear results presented in the next section. In Table 1 we also show the numerically observed convergence rate of the error measure

[TABLE]

as the time-step $h>0$ is reduced. We should note here that there seems to be a slight discrepancy in the initial data we obtain to the ones that could be seen in Fig. 7 in [13], where the same parameters are used. We obtain $|q_{1}(0)|^{2}\approx 0.1385$ and $|q_{2}(0)|^{2}\approx 0.0889$, values slightly smaller than in [13]. The reason for this is possibly different normalization.

6.5 Non-linear case

In the non-linear case we take the same initial data (6.18) as in the linear case and observe for which values of $\sigma>0$ is the beating effect supressed. All other parameters are the same as in the linear case except for

[TABLE]

We compute the solution for $\sigma=0.3,0.6,0.7,0.8,0.9,0.98$. For low values of $\sigma$ the beating phenomenon is still visible, whereas for stronger non-linearities it begins to disappear. It is also interesting that the numerical computation becomes increasingly difficult with increasing $\sigma$. This is not suprising as it is known that for large enough $\sigma$ blow-up can occur in finite time [13]. Solutions for the different $\sigma$ is given in Fugures 6–8 with some extra detail for $\sigma=0.9$ given in Figure 9. The plots for $\sigma=0.3,0.6,0.7,0.8$ do not change at this scale for smaller time-step $h$. For $\sigma=0.9$ the basic shape of the solution seems to be well captured but as indicated in Figure 9 $q_{1}$ becomes increasingly oscillatory and the error for $q_{1}$ increases significantly for larger $t$. Finally for $\sigma=0.98$ blow-up seems to occur near $t=14.3$, decreasing $h$ just increases the height of the peak. Note that the largest computation for $\sigma=0.9$ and $h=1/128$ required us to compute $N=25600$ time-steps. For this case our algorithm needed $n_{0}=57$ direct steps and $N_{Q}=112$ quadrature nodes for $\mathrm{tol}=10^{-8}$. In Figure 5 we compare the computational times of the new method, the $O(N\log^{2}N)$ method based on FFT introduced in [16] and as modified in [4, 10], and the standard naive $O(N^{2})$ implementation. These timings, clearly show that the availability of a fast method was essential to perform experiments in reasonable time. Furthermore the new method is the fastest in all the listed experiments, though the FFT based method is also very fast. However, only the new method brings savings in terms of memory. These are not significant in the one dimensional cases with a potential concentrated in only two points as investigated here. For a potential concetrated in many points, or for higher dimensional problems with the potential concentrated on a manifold, the memory savings will become equally important.

7 Conclusions

We have developed a special algorithm for the implementation of Lubich’s Convolution Quadrature when applied to the integral formulation of Schrödinger equations with concentrated potential. The new algorithm belongs to the family of the so-called fast and oblivious convolution algorithms, since for the approximation of the solution at $N$ time steps it is able to reduce the complexity from $O(N^{2})$ operations to $O(N\log N)$ and the storage from $O(N)$ to $O(n_{0}+\log N)$, with $n_{0}\ll N$. These features allow us to reliably simulate the behavior of the solution to non linear problems for long times and/or with a very small step, in order to capture high oscillations or finite time blow up. Our results are in good agreement with those reported in [21] and in [13], where two different methods with complexity $O(N^{2})$ and memory requirements $O(N)$ are used. The MATLAB codes written to perform the simulations in the current paper can be found in [5].

Future research will address the theoretical analysis of the error associated to the Convolution Quadrature approximation of relevant non linear cases, the control of the time step and the generalization of the algorithm and its application to the two dimensional case, following the recent results in [11].

Acknowledgements

The second author acknowledges Alessandro Teta for very useful discussion about the models during the preparation of the paper. The second author also acknowledges partial support by INdAM-GNCS and the Spanish grant MTM2016-75465-P.