Time adaptive Zassenhaus splittings for the Schr\"odinger equation in the semiclassical regime

TL;DR

This paper introduces a novel time-adaptive Zassenhaus splitting method for solving the Schrödinger equation in the semiclassical regime, effectively handling high oscillations and improving computational efficiency.

Contribution

It combines asymptotic Zassenhaus splitting with time adaptivity, enhancing efficiency and accuracy in numerical solutions of semiclassical Schrödinger equations.

Findings

Method reduces computational cost in high oscillation regimes

Time adaptivity optimizes step size during simulations

Numerical examples demonstrate improved accuracy and efficiency

Abstract

Time dependent Schr\"odinger equations with conservative force field U commonly constitute a major challenge in the numerical approximation, especially when they are analysed in the semiclassical regime. Extremely high oscillations originate from the semiclassical parameter, and call for appropriate methods. We propose to employ a combination of asymptotic Zassenhaus splitting with time adaptivity. While the former turns the disadvantage of the semiclassical parameter into an advantage, leading to highly efficient methods with low error constants, the latter enables to choose an optimal time step and to speed up the calculations when the oscillations subside. We support the results with numerical examples.

Figures

Tables

Table 1. Table 1: Results of different time-stepping strategies with an imposed tolerance of 10 − 7 superscript 10 7 10^{-7} for the lattice potential (left part of the table) and the Morse potential (right), ε ∈ { 10 − 2 , 10 − 3 , 10 − 4 } . 𝜀 superscript 10 2 superscript 10 3 superscript 10 4 \varepsilon\in\{10^{-2},10^{-3},10^{-4}\}. We compare the global errors with the computational effort (number of time steps respectively computation time). Adaptive strategies based on either the classical defect ‘(c)’ or the symmetrized defect ‘(s)’ are compared with fixed time steps chosen as either the smallest adaptively determined step-size, the largest adaptive step-size or the asymptotic scaling h = ε ℎ 𝜀 h=\varepsilon .

		${\psi }_L$ under lattice potential			${\psi }_M$ under Morse potential
		Global	Time	Run	Global	Time	Run
$\epsilon$	h	error	steps	time (s)	error	steps	time (s)
	adaptive (c)	$1.29\times {10}^{-6}$	490	19.5	$2.69\times {10}^{-4}$	7772	320.0
	adaptive (s)	$2.04\times {10}^{-5}$	394	18.7	$2.94\times {10}^{-4}$	7448	316.0
${10}^{-2}$	smallest (c)	$1.96\times {10}^{-7}$	746	11.9	$1.92\times {10}^{-4}$	10102	191.5
	largest (c)	$5.53\times {10}^{-4}$	217	4.1	$6.43\times {10}^{-4}$	5502	104.1
	$h=\epsilon$	$1.81\times {10}^{-2}$	100	2.4	$4.51\times {10}^{-3}$	2000	52.0
	adaptive (c)	$1.49\times {10}^{-7}$	216	12.9	$9.19\times {10}^{-4}$	13249	955.1
	adaptive (s)	$1.53\times {10}^{-7}$	216	14.3	$9.20\times {10}^{-4}$	13245	1021.5
${10}^{-3}$	smallest (c)	$1.18\times {10}^{-8}$	587	14.4	$5.00\times {10}^{-4}$	19841	572.8
	largest (c)	$1.62\times {10}^{-3}$	64	2.9	$2.21\times {10}^{-3}$	9329	315.5
	$h=\epsilon$	$3.96\times {10}^{-9}$	1000	20.4	$4.92\times {10}^{-4}$	20000	601.4
	adaptive (c)	$9.70\times {10}^{-9}$	215	43.7	$1.64\times {10}^{-3}$	27358	12029.7
	adaptive (s)	$9.75\times {10}^{-9}$	215	42.6	$1.64\times {10}^{-3}$	27343	14237.4
${10}^{-4}$	smallest (c)	$5.2\times {10}^{-10}$	925	80.2	$1.02\times {10}^{-3}$	42555	11380.8
	largest (c)	$7.12\times {10}^{-4}$	90	12.0	$3.20\times {10}^{-3}$	21307	6435.6
	$h=\epsilon$	$8.7\times {10}^{-13}$	10000	659.8	$4.95\times {10}^{-5}$	200000	60755.0

Equations

$${\partial }_t\psi (x,t)$$

$$\psi (x,0)$$

$${\psi }_1=\mathcal{S}(h){\psi }_0,$$

$$\mathcal{L}(h){\psi }_0=(\mathcal{S}(h)-\mathcal{E}(h)){\psi }_0,$$

$$\psi (h)=\mathcal{E}(h){\psi }_0$$

$$\mathcal{E}(h)={\mathrm{e}}^{\mathrm{i}h(\epsilon {\partial }_x^2-{\epsilon }^{-1}V)},$$

$$\tilde{\mathcal{L}}(h){\psi }_0=\frac{h}{p+1}\mathcal{D}(h){\psi }_0=\mathcal{O}\left(h^{p+1}\right).$$

$$h_{\mathrm{new}}=(1-\alpha )h_{\mathrm{old}}\sqrt[p+1]{\frac{\mathtt{tol}}{\tilde{\mathcal{L}}(h){\psi }_0}},\alpha >0,$$

$${\mathcal{D}}_c(t):={\partial }_t\mathcal{S}(t)-H\mathcal{S}(t)$$

$$\mathcal{L}(h){\psi }_0={\int }_0^h\mathcal{E}(h-t){\mathcal{D}}_c(t){\psi }_{0}dt=\mathcal{O}\left(h^{p+1}\right).$$

$${\tilde{\mathcal{L}}}_c(h){\psi }_0:=\frac{h}{p+1}{\mathcal{D}}_c(h){\psi }_0=\mathcal{L}(h){\psi }_0+\mathcal{O}\left(h^{p+2}\right)\text{for}h\to 0.$$

$${\partial }_t\mathcal{E}(t){\psi }_0=H\mathcal{E}(t){\psi }_0=\mathcal{E}(t)H{\psi }_0=\frac{1}{2}\{H,\mathcal{E}(t)\}{\psi }_0,$$

$${\mathcal{D}}_s(t):={\partial }_t\mathcal{S}(t)-\frac{1}{2}\{H,\mathcal{S}(t)\}.$$

$$\mathcal{L}$$

$$\mathcal{L}(h){\psi }_0={\int }_0^h\mathcal{E}(\frac{h-t}{2}){\mathcal{D}}_s(t)\mathcal{E}(\frac{h-t}{2}){\psi }_0\mathrm{d}t,$$

$${\tilde{\mathcal{L}}}_s(h){\psi }_0:=\frac{h}{p+1}{\mathcal{D}}_s(h){\psi }_0=\mathcal{L}(h){\psi }_0+\mathcal{O}\left(h^{p+2}\right).$$

$${\tilde{\mathcal{L}}}_s(h){\psi }_0=\mathcal{L}(h){\psi }_0+\mathcal{O}\left(h^{p+3}\right)\text{for}h\to 0,$$

$$\Delta x=\mathcal{O}\left(\epsilon \right),h=\mathcal{O}\left({\epsilon }^{\sigma }\right),\sigma \in (0,1].$$

$$\mathcal{S}(h)={\mathrm{e}}^{\frac{1}{2}{W}^{[0]}}{\mathrm{e}}^{\frac{1}{2}{W}^{[1]}}{\mathrm{e}}^{\frac{1}{2}{W}^{[2]}}{\mathrm{e}}^{{\mathcal{W}}^{[3]}}{\mathrm{e}}^{\frac{1}{2}{W}^{[2]}}{\mathrm{e}}^{\frac{1}{2}{W}^{[1]}}{\mathrm{e}}^{\frac{1}{2}{W}^{[0]}}=\mathcal{E}(h)+\mathcal{O}\left({\epsilon }^{7\sigma -1}\right),$$

$$W^{[0]}$$

$$W^{[1]}$$

$$W^{[2]}$$

$${\mathcal{W}}^{[3]}$$

$${⟨f⟩}_k:=\frac{1}{2}\left(f\circ {\partial }_x^k+{\partial }_x^k\circ f\right),k\ge 0(f\in {\mathrm{C}}_p^{\infty }(I;\text{mbox}\mathbb{R}))$$

$${\Vert {\mathcal{K}}_k\Vert }_2=\mathcal{O}\left((\Delta x{)}^{-k}\right)=\mathcal{O}\left({\epsilon }^{-k}\right).$$

$$\widetilde{{⟨f⟩}_k}=\frac{1}{2}\left({\mathcal{D}}_f{\mathcal{K}}^k+{\mathcal{K}}^k{\mathcal{D}}_f\right),$$

$${\mathcal{S}}_j(t)=\exp (t^{n_j}{\mathcal{R}}^{[j]}),$$

$${\partial }_t{\mathcal{S}}_j(t)=n_j{t}^{n_j-1}{\mathcal{R}}^{[j]}\exp (t^{n_j}{\mathcal{R}}^{[j]}).$$

$$\mathcal{S}(t)$$

$${\mathcal{R}}^{[0]}$$

$${\mathcal{R}}^{[1]}$$

$${\mathcal{R}}^{[2]}$$

$${\mathcal{R}}^{[3]}$$

$${\partial }_t\mathcal{S}(t)$$

$$v_1$$

$$v_2$$

$$v_3$$

$$({\psi }_1=)v_4$$

Full text

††journal: Appl. Math. Comput.

url]http://www.asc.tuwien.ac.at/ winfried

url]http://www.harald-hofstaetter.at

url]https://mat.ug.edu.pl/ kmalina

url]https://www.maths.ox.ac.uk/people/pranav.singh

1 Introduction

In this paper we are concerned with developing a time-adaptive method for solving the Schrödinger equation in the semiclassical regime,

[TABLE]

with periodic boundary conditions imposed on the interval $I\subseteq\mbox{\Bbb R}$. The interaction potential $V(x)$ is a real and periodic function. The regularity required of $V$ and $\psi$ depends on the desired order of the numerical method. For the sake of simplicity, we assume $V\in\mathrm{C}_{\mathrm{p}}^{\infty}(I;\mbox{\Bbb R})$ and $\psi\in\mathrm{C}_{\mathrm{p}}^{\infty}(I;\mbox{\Bbb C})$, where the subscript “$\mathrm{p}$” denotes periodicity.

The semiclassical parameter $\varepsilon$ induces oscillations of wavelength ${\mathcal{O}}\!\left(\varepsilon\right)$ in the solution $\psi$, both in space and in time [1, 2, 3]. In this regime, finite difference schemes require an excessively fine spatial grid and very small time steps [4]. Consequently, they are found to be ineffective in comparison with spectral discretisation in space followed by exponential splittings for time-propagation [1]. Methods based on Lanczos iterations [5], which are particularly effective for problems in the atomic scaling ($\varepsilon=1$), are also found to be ineffective in the semiclassical regime, where the exponent involved is of huge spectral size, scaling as ${\mathcal{O}}\!\left(\varepsilon^{-1}\right)$, see [6, 7]. Effective methods for highly oscillatory problems, which have a Hamiltonian structure and are periodic in time, have been proposed in [8] and are referred to as Hamiltonian boundary value methods. They have been applied in the context of Schrödinger equations in [9].

In this regime the symmetric Zassenhaus splittings of [10] are found to be very effective. These are asymptotic exponential splittings where exponents scale with powers of the small parameter $\varepsilon$ and, consequently, become progressively small. The small size of these exponents allows a very effective approximation via Lanczos iterations despite reasonably large time steps [10].

Needless to say, the oscillations in the solution change in time. Higher oscillations require smaller time steps, which increases the computational cost. Thus, we face the usual tradeoff between two competing concerns – smaller time steps for higher accuracy and larger time steps for lower cost. Time adaptivity has been a well developed approach to arrive at the optimal compromise by keeping the time steps as large as possible, so long as a prescribed error tolerance is not exceeded. State-of-the-art step-size choice is based on firm theoretical ground by recent investigations of digital filters from signal processing and control theory [11, 12, 13], which have also been demonstrated to enhance computational stability [14]. The advantages of adaptive selection of the time steps in the context of splitting methods have been demonstrated for nonlinear Schrödinger equations in [15], and for parabolic equations in [16]. It is found that in addition to a potential increase in the computational efficiency, the reliability of the computations is enhanced. In an adaptive procedure, the step-sizes are not guessed a priori but chosen automatically as mandated by the smoothness of the solution, and the solution accuracy can be guaranteed.

The aim of this paper is thus to design a time-adaptive method for the Schrödinger equation in the semiclassical regime. This is done by utilising the defect based time adaptivity approach of [17, 18] for the symmetric Zassenhaus splittings of [10]. These time adaptivity schemes, which are very effective for classical splittings, turn out to be successful for the asymptotic splittings of [10] as well.

In Section 2 we describe the defect-based time-adaptive approach. In Section 3 we present a variant of the high-order Zassenhaus splittings of [10], which turns out to be more conducive in a time-adaptive approach. The practical algorithm used for the estimation of the local error for this highly efficient method is derived in Section 4. In Section 5 we present numerical experiments which demonstrate the efficacy of the time-adaptive scheme.

2 Defect-based time adaptivity

Time-adaptivity in numerical solutions of ODEs and PDEs involves adjusting time step of a numerical scheme in order to keep the local error in a single step below a specified error tolerance, ${\tt tol}$. The local error in a single step of a one-step integrator with step-size $h$, starting (without loss of generality) from $\psi_{0}$ at $t=0$,

[TABLE]

where $\mathcal{S}(h)$ is the numerical evolution operator, is given by

[TABLE]

where

[TABLE]

is the exact solution. Since the exact evolution operator for (1.1),

[TABLE]

and the exact solution are not available, practical time-adaptivity methods rely on accurate a posteriori estimates of the local error, $\tilde{\mathcal{L}}(h)\psi_{0}$, that can be computed along with the numerical solution. In this manuscript, we will focus on defect-based estimators of the local error, which are of the form

[TABLE]

Here, $\mathcal{D}(h)\psi_{0}$ is a computable defect term, measuring the local quality of the approximation delivered by the numerical solution, $\mathcal{S}(h)\psi_{0}$. We will consider two different versions of the defect $\mathcal{D}(h)\psi_{0}$, see Sections 2.1 and 2.2.

Once an estimate of the local error via (2.1) is available, a new step-size can then be chosen as

[TABLE]

where ${\tt tol}$ is the local error tolerance and $p+1$ is the local order of the numerical scheme, $\mathcal{S}(h)=\mathcal{E}(h)+{\mathcal{O}}\!\left(h^{p+1}\right)$. The factor of $(1-\alpha)$ ensures that we are more conservative with the time step, always taking a slightly smaller time step than predicted ($\alpha$ is usually taken to be a small positive number such as $0.1$). When the local error estimate in a step $\widetilde{\mathcal{L}}(h)\psi_{0}$ exceeds the error tolerance (${\tt tol}$), the numerical propagation is run again with the smaller step $h=h_{\mathrm{new}}$, otherwise the new time step is used for the next step.

2.1 Classical defect-based estimator

We briefly recall the idea underlying (2.1); see for instance [17]. Thinking of the step-size $h$ as a continuous variable and denoting it by $t$ (somewhat more intuitive), the discrete flow $\mathcal{S}(t)\psi_{0}$ is a well-defined, smooth function of $t$. We call

[TABLE]

the classical defect (operator) associated with $\mathcal{S}$, obtained by plugging in the numerical flow into the given evolution equation (1.1) (which is satisfied exactly by its exact flow $\mathcal{E}(t)$). Due to the definition of $\mathcal{D}_{c}$, the local error enjoys the integral representation

[TABLE]

Evaluation of ${\widetilde{\mathcal{L}}}_{c}(h)\psi_{0}$ requires a single evaluation of the defect at $t=h$, the step-size actually used in the computation of $\psi_{1}=\mathcal{S}(h)\psi_{0}$.222How to compute the defect at $t=h$, for the given step-size $h$ used in the computation, in the context of Zassenhaus splitting will be explained in Section 4.

2.2 Symmetrized defect-based estimator

An alternative way of defining the defect was introduced in [18, 19]. It is based on the fact that the exact evolution operator $\mathcal{E}(t)$ commutes with the Hamiltionian $H$, whence

[TABLE]

with the anti-commutator $\{H,X\}=HX+XH$. This motivates the definition of the symmetrized defect (operator)

[TABLE]

[TABLE]

Now suppose that $\mathcal{S}$ is symmetric (time-reversible), i.e., it satisfies ${\mathcal{S}(-t)\mathcal{S}(t)}=I$. Then, its order $p$ is necessarily even, and moreover we even have

[TABLE]

see [18, 19]. This means that in the symmetric case the deviation $({\widetilde{\mathcal{L}}}_{s}(h)-\mathcal{L}(h))\psi_{0}$ of the local error estimate based on the symmetrized defect is of a better quality, asymptotically for $h\to 0$, than that one based on the classical defect.333Note that $\mathcal{D}_{c}(h)={\mathcal{O}}\!\left(h^{p}\right)$, $\mathcal{D}_{s}(h)={\mathcal{O}}\!\left(h^{p}\right)$, and $\mathcal{D}_{s}(h)-\mathcal{D}_{c}(h)={\mathcal{O}}\!\left(h^{p+1}\right)$. Moreover, evaluation of $\mathcal{D}_{s}(h)\psi_{0}$ is typically only slightly more expensive compared to $\mathcal{D}_{c}(h)\psi_{0}$, see [18, 19] for the case of conventional splittings and, in particular, Section 4 below.

3 Symmetric Zassenhaus splittings

Symmetric Zassenhaus splittings are asymptotic splittings introduced in [10] for solving the Schrödinger equation in the semiclassical regime (1.1). These splittings are derived by working in infinite dimensional space, prior to spatial discretisation, which enables a more effective exploitation of commutators, since some of them vanish, while the rest end up being much smaller (in a sense of spectral radius) than generically expected.

To manage the ${\mathcal{O}}\!\left(\varepsilon\right)$ wavelength oscillations in space and time, which arise due to the presence of the small semiclassical parameter $\varepsilon$, the following relations have been established to be useful for the spatio-temporal resolution in these schemes.

[TABLE]

The Zassenhaus splitting that we will consider in this paper is

[TABLE]

Here,

[TABLE]

are the symmetrized differential operators which first appeared in [10] and have been studied in detail in [20, 3]. These operators are used for preservation of stability under discretization after simplification of commutators.

Once a differential operator such as $\partial_{x}^{2}$ is discretised to a symmetric differentiation matrix $\mathcal{K}_{2}$ via spectral collocation, its spectral size grows as $\Delta x={\mathcal{O}}\!\left(\varepsilon\right)$ decreases, since

[TABLE]

Keeping eventual discretisation with the scaling (3.1) in mind, we use a shorthand $\partial_{x}^{k}={\mathcal{O}}\!\left(\varepsilon^{-k}\right)$ for the undiscretised and unbounded operator as well. The symmetrized differential operator $\left\langle f\right\rangle_{k}$ discretises to the form

[TABLE]

where $\mathcal{D}_{f}$ is a diagonal matrix with values of $f$ at the grid points. Its spectral radius grows as ${\mathcal{O}}\!\left(\varepsilon^{-k}\right)$, assuming that $f$ is independent of $\varepsilon$, and we abuse notations once more to say $\left\langle f\right\rangle_{k}={\mathcal{O}}\!\left(\varepsilon^{-k}\right)$, for short.

The splitting (3.2a) possesses several favorable features. First of all, the critical quantities like time step $h$, spatial step $\Delta x$ and semiclassical parameter $\varepsilon$ are tied together by relation (3.1). As a consequence, the splitting error is expressed via a universal quantity $\mathcal{O}(\varepsilon^{7\sigma-1})$, and the error constant does not hide any critical quantities.

The exponentials involved in the splitting are easily computable after spatial discretization. For example, spectral collocation transforms $W^{[0]}$ into a symmetric circulant matrix, whence its exponential may be computed by the Fast Fourier Transform. The diagonal matrix $W^{[1]}$ can be exponentiated directly. Neither $W^{[2]}$ nor $\mathcal{W}^{[3]}$ are structured favourably, however their spectral radius is small enough to exponentiate them with a small number (say, 3 or 4) of Lanczos iterations.

Remark 1

Since $h={\mathcal{O}}\!\left(\varepsilon^{\sigma}\right)$ the splitting (3.2a) which features an error of ${\mathcal{O}}\!\left(\varepsilon^{7\sigma-1}\right)$ could loosely be considered a sixth-order splitting with an error constant scaling as $\varepsilon^{-1}$. This interpretation is not strictly correct, however, since the error estimate holds for $\sigma\in(0,1]$ and $\varepsilon\rightarrow 0$. For a finite $\varepsilon=\varepsilon_{0}$ and $h\rightarrow 0$ (i.e. $\sigma\rightarrow\infty$), the asymptotic Zassenhaus splitting (3.2a) is a highly efficient fourth-order method with a very small error constant. This disparity in the different asymptotic limits occurs because the derivation of (3.2a) involves discarding terms of size ${\mathcal{O}}\!\left(h^{5}\varepsilon\right)$, for instance. These terms are ${\mathcal{O}}\!\left(\varepsilon^{7\sigma-1}\right)$ under the asymptotic scaling (3.1), $\varepsilon\rightarrow 0,h={\mathcal{O}}\!\left(\varepsilon^{\sigma}\right),\sigma\leq 1$, but are ${\mathcal{O}}\!\left(h^{5}\right)$ under a fixed $\varepsilon=\varepsilon_{0}$ and $h\rightarrow 0$. The consequence for the time-adaptivity algorithm (2.1) is that we must use $p=4$, not $p=6$.

Remark 2

Many different variants of this splitting are possible (including, but not limited to the case when we start with $W^{[0]}=-{\mathrm{i}}\,h\,\varepsilon^{-1}V$). Application of these variants can just as easily be considered, however will not be the focus of the present paper.

Remark 3

The exponents (3.2b) differ from the splitting described in [10] in one minor aspect – the term $-\tfrac{1}{24}{\mathrm{i}}\,h^{3}\varepsilon(\partial_{x}^{4}V)$ has been moved from $\mathcal{W}^{[3]}$ to $W^{[2]}$. This term is of size ${\mathcal{O}}\!\left(\varepsilon^{3\sigma+1}\right)$, which is smaller than ${\mathcal{O}}\!\left(\varepsilon^{3\sigma-1}\right)$ and thus may also be combined with $W^{[2]}={\mathcal{O}}\!\left(\varepsilon^{3\sigma-1}\right)$. In [10] this term is combined with $\mathcal{W}^{[3]}={\mathcal{O}}\!\left(\varepsilon^{5\sigma-1}\right)$ working under the assumption $\sigma\leq 1$ (time steps larger than ${\mathcal{O}}\!\left(\varepsilon\right)$) since ${\mathcal{O}}\!\left(\varepsilon^{3\sigma+1}\right)$ is also smaller than ${\mathcal{O}}\!\left(\varepsilon^{5\sigma-1}\right)$ under $\sigma\leq 1$. Consequently, all exponents feature a single power of $h$. This change makes no difference to the overall order of the scheme but makes the computation of defect easier.

4 Local error estimator for Zassenhaus splitting

Remark 4

In the sequel, the argument $(t)$ or $(h)$, respectively, is suppressed whenever the meaning is obvious.

4.1 Time derivative of the discrete evolution operator

The time-adaptivity algorithm is based on the estimation of the local error in each step, which in turn requires the computation of the defect at discrete time $t=h$, see Section 2. This requires the computation of $\partial_{t}\mathcal{S}(h)=\partial_{t}\mathcal{S}(t)|_{t=h}$ For splitting schemes such as (3.2a), this boils down to an application of the product rule, involving computation of the time derivative of each individual exponential.

For exponentials such as

[TABLE]

Now, consider the (modified) Zassenhaus splitting (3.2a) & (3.2b) for (1.1), with the discrete evolution operator re-written in the form

[TABLE]

where $\mathcal{S}_{j}=\mathcal{S}_{j}(t)$, $n_{0}=1,\ n_{1}=1,\ n_{2}=3,\ n_{3}=5$, and

[TABLE]

The derivative of the flow for a splitting scheme can be expressed via the product rule,

[TABLE]

Here the underlined exponentials are the ones that need to be computed freshly. The rest can be dealt with by storing intermediate values from evaluation of the splitting scheme: In order to compute $\partial_{t}\mathcal{S}(h)$, we store

[TABLE]

which are anyway computed during evaluation of the numerical scheme. In addition, we compute and store

[TABLE]

Then,

[TABLE]

Thus, we need to compute six exponentials appearing in (4.5) in addition to the seven exponentials required in the Zassenhaus splitting (3.2a). For a scheme of order six, we need a total of 13 exponentials.

Remark 5

For the order four method

[TABLE]

five exponentials are required for the Zassenhaus splitting and we have verified that four additional exponentials are required for time adaptivity, making for a total of nine exponentials.

Remark 6

In a practical, memory-efficient implementation, the $w_{j}$ are computed ‘on the fly’ along with the $v_{j}$, via alternating updates of arrays $v$ and $w$.

Remark 7

In the above computations, the only feature specific to the Zassenhaus splitting is the use of (4.1b) for the derivative of the exponential. Thus the derivation and the combination of exponentials works in an analogous way for other splittings, with different time derivatives.

4.2 Practical evaluation of the classical defect (2.2)

We need to compute the classical defect at discrete time $t=h$,

[TABLE]

Here, $\psi_{1}=\mathcal{S}(h)\psi_{0}$, $H=2(\mathcal{R}^{[0]}+\mathcal{R}^{[1]})$ is the Hamiltonian, and $\partial_{t}\mathcal{S}(h)\psi_{0}$ is the time derivative according to Section 4.1, see (4.5).

In addition, a small amount of work can be saved by exploiting the relation $w_{4}=\mathcal{R}^{[0]}v_{4}$ (see 4.4). Then, computation of $\mathcal{D}_{c}(h)\psi_{0}$ amounts to the evaluation of

[TABLE]

4.3 Practical evaluation of the symmetrized defect (2.4)

Instead of (4.6), we need to compute

[TABLE]

With (4.5), we have (again making use of $w_{4}=\mathcal{R}^{[0]}v_{4}$)

[TABLE]

Here the underlined exponential is the one that needs to be computed in addition, compared to evaluation of $\mathcal{D}_{c}(h)\psi_{0}$. Of course, this additional effort is marginal.

Thus, computation of the symmetrized defect at discrete time $t=h$ amounts to evaluation of

[TABLE]

5 Numerical experiments

For our numerical experiments, we will use the wave-packets

[TABLE]

as initial conditions for two experiments. Here

[TABLE]

is a wave-packet with a spread of $\delta$, mean position $x_{0}$ and mean momentum $k_{0}$. We will consider the evolution of $\psi_{L}$ as it heads towards the lattice potential $V_{L}$ and the evolution of $\psi_{M}$ as it oscillates in the Morse potential $V_{M}$, which are given by

[TABLE]

respectively, where

[TABLE]

is a bump function.

We consider the behaviour at the final time $T_{L}=1$ for the first experiment and $T_{M}=20$ for the second experiment under the choices of $\varepsilon=10^{-2},10^{-3}$ and $10^{-4}$. The spatial domain is chosen as $[-2,2]$ for the lattice potential and $[3,10]$ for the Morse potential, and we impose periodic boundary conditions. The evolution of these wavepackets under the choice $\varepsilon=10^{-2}$ is shown in Fig. 1.

As one can observe in Figures 2 and 3, the procedure quickly adapts the time step to ensure that the local accuracy is within the specified threshold, which is taken to be ${\tt tol}=10^{-7}$ here.

In the first experiment, we use $M=750,1750$ and $15000$ equispaced points for spatial discretization in the cases $\varepsilon=10^{-2},10^{-3}$ and $10^{-4}$, respectively, while $M=500,1500$ and $10000$ are used for the three choices of $\varepsilon$ in the second experiment.

In Table 1, we show the total number of time steps required for maintaining the local accuracy of $10^{-7}$ for three cases: $\varepsilon=10^{-2},10^{-3}$ and $10^{-4}$. Also presented are the global $\mathrm{L}^{2}$ accuracies of these solutions at the final time ($T_{L}=1$ and $T_{M}=20$) and the computational time. These are compared to the (non-adaptive) case when the time step is fixed at (i) the finest value used by the adaptive method (ii) the coarsest value used by the adaptive method (iii) $h=\varepsilon$. The reference solutions for computing global errors are produced by using Matlab’s expm function, while using a much finer spatial grid. The number of Lanczos iterations are chosen automatically to meet the accuracy requirements using the a priori error bounds of [21].

We find that for small $\varepsilon$, the asymptotic scaling $h=\varepsilon$ yields the smallest error, but the computational effort is prohibitive. The symmetrized error estimate invariantly generates coarser time grids than the classical version, but the additional computational effort implies that the picture is ambivalent in that the computational effort is smaller in only about half of the cases. When the largest step generated by the adaptive strategy ‘(c)’ is used in a uniform grid, computations are fast but not sufficiently accurate, while use of the smallest adaptive time step is prohibitively expensive in the case of the lattice potential (for the Morse potential, the reduction in the number of time steps is outweighed by the computational effort for the error estimate). Together, these two observations imply that adaptive step-size choice constitutes the appropriate means to determine grids that reproduce the solution in a reliable and efficient way, and excels over a priori choice of the step-size. In particular, the guess of an optimal step-size is commonly not feasible, while the adaptive strategy finds the time steps automatically and guarantees the desired level of accuracy.

6 Conclusions

We have investigated adaptive strategies used in conjunction with asymptotic Zassenhaus splitting for the solution of linear time-dependent Schrödinger equations in the semiclassical regime. The local time-stepping is based on defect-based error estimation strategies. We have demonstrated that adaptivity provides a means to reliably and efficiently achieve a desired level of accuracy especially for a small semiclassical parameter and thus excels over fixed time steps in many cases. A major benefit is in addition that no near-optimal equidistant step-size has to be guessed a priori.