Fourier methods for the perturbed harmonic oscillator in linear and nonlinear Schrödinger equations

0.29 MB·English

by  Philipp Bader

Fourier methods for the perturbed harmonic oscillator in linear and nonlinear
Schr¨odinger equations
Philipp Bader and Sergio Blanes
∗ †
Universitat Polit`ecnica de Val`encia, Instituto de Matema´tica Multidisciplinar, E-46022 Valencia, Spain
WeconsiderthenumericalintegrationoftheGross-Pitaevskiiequationwithapotentialtrapgiven
by a time-dependent harmonic potential or a small perturbation thereof. Splitting methods are
frequentlyusedwithFouriertechniquessincethesystemcanbesplitintothekineticandremaining
part,andeachpartcanbesolvedefficientlyusingFastFourierTransforms. Tosplitthesysteminto
thequantumharmonicoscillator problemandtheremainingpartallows togethigheraccuracies in
many cases, but it requires to change between Hermite basis functions and the coordinate space,
1 and this is not efficient for time-dependent frequencies or strong nonlinearities. We show how to
1 buildnewmethodswhich combinetheadvantagesofusingFouriermethodswhilesolving thetime-
0 dependent harmonic oscillator exactly (or with a high accuracy by using a Magnus integrator and
2 an appropriate decomposition).

n
a PACSnumbers: 02.60.-x,02.60.Cb,02.60.Lj,02.70.Hm
J

1
3 I. INTRODUCTION accuracy with a moderate number of mesh points and
low computational cost. For harmonic oscillator (HO)

] The numerical integration of the Gross-Pitaevskii problems, however, the exact solution is known and by
A
equation (GPE) expanding the solution in Hermite polynomials, highly
N accurate results are obtained if the HO is solved sepa-
h. i ∂ ψ(x,t)= 1 ∆+V(x,t)+σ(t)ψ(x,t)2 ψ(x,t) rately [2, 3, 6, 7].
∂t −2µ | | It is claimed [2], that either Hermite or Fourier pseu-
t (cid:18) (cid:19)
a dospectralmethods are the most efficient, the choice de-
m x Rd,describingthegroundstateofinteractingbosons
pending on the particular parameter set. Motivated by
∈
at zero temperature, the Bose-Einstein condensates, has
[ theseresults,weshowhowbothmethodsarecombinedto
attracted great interest [1–3] after the first experimental
retainboth the accuracyofthe Hermite method and the
2 realizations[4]. Wepresentanewefficientwaytosolvea
speed of the Fourier transforms, i.e. to rewrite the Her-
v
special class of GPE, namely that of weakly interacting
0 mite method as a single simple pseudospectral Fourier
bosons in a single time-dependent trap. To be more spe-
7 scheme. We have found, that this approximation per-
cific,thepotentialtrapV istakentobeaperturbationof
4 forms, for the studied problem class, always equal to or
3 the (time-dependent) d-dimensional harmonic oscillator, better than the original Fourier method and therefore
7. i.e. V(x,t) = xTM(t)x+ǫVI(x,t) where M(t) ∈ Rd×d has to compete with Hermite expansions only. Hermite
is a positive definite matrix and ǫV (x,t) is a small per-
0 I schemes suffer from large computational costs when the
0 turbation. Therealscalarfunctionσ originatesfromthe number of basis terms in the expansions is altered along
1 mean-field interaction between the particles and corre- the integration or taken very large, being the case for
: sponds to repulsive or attractive forces for positive or
v time dependent trap frequencies M(t) or strong nonlin-
negative values of σ(t), respectively [5]. Notice that the
i earities σ(t). It is in this setup, where our new method
X non-interacting case, σ 0, corresponds to the linear
substantially improves the Hermite performance, and it
≡
r Schr¨odinger equation. canindeedberegardedastheoptimalchoiceforthenum-
a
Several methods have been analyzed to compute both
ber of Hermite basis functions (for an equidistant grid)
the time evolution and the ground state of the GPE in
at each time step.
thecourseofthelastdecade[1–3,6,7],amongthemfinite
For the ease of notation, we restrict ourselves to the
differences, Galerkin spectral methods and pseudospec-
one-dimensional problem
tral methods for Fourier or Hermite basis expansions. It
has been concluded [2] that these pseudospectral meth- ∂
i ψ =H (t)ψ+ εV (x,t)+σ(t)ψ 2 ψ (I.1)
odsperformbestforawideparameterrangefortheGPE. ∂t 0 I | |
TheFouriertypemethodscanbeimplementedwithFast
(cid:0) (cid:1)
Fourier Transform (FFT) algorithms since the trapping where
V causesthewavefunctiontovanishasymptoticallythus
1 1
allowing to consider the problem as periodic on a suffi- H0(t)= p2+ µω2(t)x2 (I.2)
2µ 2
ciently large spatial interval. Their advantages are high
andp= i ∂ . The boundaryconditionsimposedbythe
− ∂x
trap require the wave function to go to zero at infinity,
∗ [email protected] and up to any desired accuracy, we can assume ψ(x,t)
† [email protected] andallitsderivativestovanishoutsideafiniteregion,say
2
[a,b], which we divide using a mesh (usually with N = thateliminated one computationis calledFirstSame As
2k points to allow a simple use of the FFT algorithms). Last (FSAL) property and can also be used with higher
Then, the partial differential equation (I.1) transforms order m stage compositions
−
into a system of ordinary differential equations (ODEs)
ψ[p] =ϕ[A] ϕ[B] ϕ[A] ϕ[B] (II.5)
d h amh◦ bmh◦···◦ a1h◦ b1h
i u(t)=Hˆ (t)u(t)+ εV (X,t)+σ(t)u(t)2 u(t),
dt 0 I | | if am = 0 or bm = 0 and repeated application of the
(cid:0) (cid:1) (I.3) scheme without requiring output. For linear problems,
and the harmonic part becomes it is usual to replace the flow-maps by exponentials (for
nonlinear problems the same is possible using exponen-
1
Hˆ (t)=T + µω2(t)X2, (I.4) tials of Lie operators). In this notation, the equation
0
2 iu = A(u)+B(u), whose formal solution for the evolu-
′
with u CN, where ui(t) ψ(xi,t), xi = a+ih, i = tion operator is denoted by φt[A+B] = e−it(A+B), is ap-
0,1,...,∈N 1, h=(b a)≃/N, X =diag x ,...,x proximated for one time step, h, by the order p compo-
0 N 1
and T deno−tes a discret−ization of the kine{tic part. − } sition (II.5) or, equivalently,
This system of nonlinear ODEs can be numerically
solved by standard all purpose ODE methods. However, ψh[p] ≡e−ihamAe−ihbmB ··· e−iha1Ae−ihb1B. (II.6)
because of the particular structure of this problem, dif-
We keep in mind that, in a nonlinear problem, if B
ferent numerical methods can differ considerably in ac-
dependsonu,ithastobeupdatedateachstagebecause
curacy as well as in computational cost and stability. In
u changes during the evolution of e ihaiA.
addition, the structural properties of the system lead to −
Backward error analysis shows that the action of a
theexistenceofseveralpreservedquantitieslikethenorm
splitting method is equivalent to solve exactly, for one
and energy (for the autonomous case).
time step h, a perturbed differential equation (for suffi-
Theaccuratepreservationofthesequantitiesaswellas
ciently small h) that can contain higher order commuta-
theerrorpropagationandperformanceofsplittingmeth-
tors
odsexplainwhytheyarefrequentlyrecommendedforthe
time integration [2, 3, 6] and make them subject of in-
u =(A+B+hα [A,B] (II.7)
′ 1,1
vestigation in this work.
+h2(α [A,[A,B]]+α [B,[B,A]])+... u
2,1 2,2
([A,B]:=AB BA,andforsimplicity, wehaved(cid:1)enoted
II. SPLITTING METHODS −
A=A(u) , B =B(u) ) where the coefficients a ,b
i i
·∇ ·∇
are chosen in order to cancel out the coefficients α up
i,j
Let us consider the separable system of ODEs
toagivenorder. Itisthusimportanttoanalyzethedom-
u =A(u)+B(u), u(t )=u CN, (II.1) inant contributions to the error from the commutators,
′ 0 0 ∈ in order to choose the most appropriate method or to
where we assume that both systems build new ones.
Thereexistmanydifferentsplittingmethodswhichare
u′ =A(u), u′ =B(u) (II.2) designed for different purposes, depending on the struc-
ture of the problem, the desired order, the required sta-
can either be solved in closed form or accurately inte- bility, etc. [8–14]. If the operator A corresponds to the
grated. If ϕ[A], ϕ[B] represent the exact flows associated kinetic part,whichis quadraticinmomenta,andthe op-
t t
to (II.2) then, to advance the solution one time step, h, erator B is the potential, diagonal in coordinate space,
wecanuse,forexample,thecompositionψ[1] =ϕ[A] ϕ[B] the commutator [B,[B,[B,A]]] vanishes and partitioned
h h ◦ h Runge-Kutta-Nystr¨om methods become favorable. The
[1] [A] [B]
(i.e. u(t +h) ψ (u ) = ϕ ϕ (u ) ), which is
0 ≃ h 0 h h 0 coefficients for this family of methods have to solve a
known as the first-order Lie-Trotte(cid:16)r method(cid:17). A method significantly reduced number of order conditions and, in
has order p if ψh[p] = ϕh + O(hp+1) where ϕt denotes general, their performance is superior to splitting meth-
the exact global flow of (II.1). Sequential application of ods designed for general separable problems, in addition
the two first-order methods ψ[1] and its adjoint ψ[1]∗ = to being more stable. Efficient schemes of order 4 and 6
h h
ϕ[B] ϕ[A] with half time step yields the second order are obtained in [9]. On the other hand, when H0 is the
h ◦ h dominant part, it is worth to take a closer look at the
time-symmetric methods
split(I.1). Thiswouldcorrespondto B A andfor
k k≪k k
ψ[2] =ϕ[A] ϕ[B] ϕ[A] (II.3) this case, when facing autonomous problems, there exist
h,A h/2◦ h ◦ h/2 tailored methods which have shown a high performance
ψ[2] =ϕ[B] ϕ[A] ϕ[B] (II.4) in practice. Writing the equation as iu = (A +ǫB)u
h,B h/2◦ h ◦ h/2 ′
(with ǫ a small parameter), it is clear that the local er-
(referredasABAandBAB compositions). Thecontrac- ror of the second order methods (II.3) or (II.4) comes
tion via the (1-parameter-)group property of the flows from the commutators at third order ([A,[A,ǫB]] and
3
[ǫB,[A,ǫB]]) and we can say that the local error is of H of (I.1),
0
order (ǫh3+ǫ2h3). Thecoefficientsa ,b inthe general
i i
O
composition (II.6) can be chosento cancel the dominant ∂ 1 1
i ψ = p2+ µω2(t)x2 ψ (II.11)
errorterms,say,the (ǫhr)termsforrelativelylargeval- ∂t 2µ 2
O (cid:18) (cid:19)
ues of r. Then, one can denote the effective order of a
method by (r,p) with r p when the localerroris given is easily computed for a time step using Fourier trans-
≥
by (ǫhr+1 +ǫ2hp+1). The method is of order p, but forms. Before giving the details on the time integration,
O
in the limit ǫ 0 it is considered as of order r p. some remarks on the formal solution are necessary.
→ ≥
Using this split allows to gain a factor ǫ in the accuracy It is well known that H is an element of the Lie alge-
0
even for general splitting methods where r = p. In [11], bra spanned by the operators E =x2/2,F =p2/2,G=
several methods of order (r,2) for r 10 are obtained 1(px+xp) , where µ=1 for s{implicity, and its commu-
with all coefficients ai,bi being positiv≤e and some other t2ators are }
schemes of order (r,4) for r = 6,8 are presented. For
near-integrablesystems,theselastmethods arethe most [E,F]=iG, [E,G]=2iE, [F,G]= 2iF.
efficient and stable ones. Despite the gain of accuracy, −
the split into a dominant part and a small perturbation This is a three-dimensionalLie algebraandthe solution,
is left unconsidered when it leads to involved or compu- ψ(x,t) = U(t,0)ψ(x,0), of (II.11) can be expressed as a
tationally costly algorithms. This issue is addressed in singleexponentialusingtheMagnusseriesexpansion[16,
this work. 17] or as a product of exponentials [18]. It is possible to
To take accountfor the time dependence in the vector formulatetheevolutionoperatorU(t,0)inmanydifferent
fields in (I.1), a more detailed analysis is required. The
ways, the most appropriate depending on the particular
general separable equation purpose, e.g. using the Magnus expansion
iu =A(u,t)+ǫB(u,t) (II.8)
′ U(t,0)=exp(f (t)E+f (t)F +f (t)G) (II.12)
1 2 3
can be solved by considering the time as two new inde-
forcertainfunctionsf (t)[19]. Approximationsof (II.12)
pendent coordinates i
for one time step, h, on the other hand, are easily ob-
iu = A(u,t ) iu = ǫB(u,t ) tained, e.g. a fourth-order commutator-free method is
′ 2 , ′ 1 . (II.9)
t = 1 t = 1 given by [20]
(cid:26) ′1 (cid:26) ′2
Thisstandardsplitisequivalenttoanewsysteminan h 1 1
U(t+h,t)=exp ( p2+ω2 x2) (II.13)
extended phase space which is not near-integrable and 2 2 L2
thehighlyefficientandstablemethodsdesignedforthose (cid:18) (cid:19)
h 1 1
problemslosetheirexcellentperformance. Following[15], exp ( p2+ω2 x2) + (h5)
× 2 2 R2 O
thenearintegrabilityisrecoveredifweintroduceonlyone (cid:18) (cid:19)
variable as follows
where
iu = A(u,t )
(cid:26) t′1′ = 1 1 , iu′ =ǫB(u,t1). (II.10) ωL2 =αω12+βω22, ωR2 =βω12+αω22
This split requires to solve exactly the equation iu′ = with ω = ω(t + c h), c = 1 √3, c = 1 + √3,
A(u,t)forthecorrespondingfractionaltimestepsandto i n i 1 2 − 6 2 2 6
and α = 1 1 , β = 1 α. It can be considered
freezethetimewhensolvingiu′ =ǫB(u,t1)(seeTableI). 2 − √3 −
as the composition of the evolution for half time step
of two oscillators with averaged frequencies, using the
TABLEI.Algorithmforthenumericalintegrationofthesys- fourth-orderGauss-Legendrequadratureruletoevaluate
tem (II.10) bythe composition (II.5). ω(t). Different quadrature rules can also be be used and
Algorithm for one time step t →t +h correspond to different averages along the time step, see
n n
[17,20]. Inthelimitwhenωisconstanttheexactsolution
t[0]=t
n
is recovered. Higher order approximations are available,
do i=1,m
ifmoreaccurateresultsaredesired,byapproximatingthe
solve: u′ =B(u,t[i−1]), t∈[t[i−1],t[i−1]+b h]
i functions f in (II.12) via truncated Magnus expansions.
i
solve: u′ =A(u,t), t∈[t[i−1],t[i−1]+a h]
i Our objective is to obtain a factorization of the solu-
t[i] =t[i−1]+aih tion which only involves terms proportional to E or F
enddo since they are easy to compute as we show in the follow-
t =t[m] ing paragraph.
n+1
Starting from (II.13) or high order approximates of
We show that the exact solution of the non- (II.12),themainresultofthisworkisthenaturaldecom-
autonomous problem, in our setting the dominant part position for the application of Fourier spectral methods.
4
A. Hamiltonians for spectral methods The previous results show how to compute individual
parts of the equation and thus permit different ways of
We now analyze how to compute the evolution of dif- splitting the system in two solvable parts
ferent parts of the Hamiltonian by spectral methods.
iψ =(A+B)ψ (II.20)
The spatial derivative (or kinetic part) associatedto the t
semidiscretized problem (I.3) can be solved in the mo- We consider the following cases:
mentum space by noting that (i) Fourier(F)-split.
d
idtu(t)=Tu(t)=FN−1DNFNu, (II.14) A= 1p2, B(t)= ω(t)x2+εV (x,t)+σ(t)ψ 2,
I
2 2 | |
where DN is diagonal and N denotes the discrete (II.21)
F
Fourier transform of length N, whose computation can We take the time as a new coordinate, as in (II.10), and
be accomplishedby the FFT algorithmwith (NlogN) evolve it with A, which is now autonomous and exactly
O
floating point operations. The solution of (II.14) for one solvable, and freeze the time in B, which is then solved
time step, h, is given by using the result from Lemma II.1. Here, A and B are
diagonal in the momentum and coordinate spaces, re-
u(t+h)=FN−1e−ihDNFNu(t) spectively, and we can change between them using the
Fourier Transforms.
which requires two FFT calls. The exponentials in
(ii) Harmonic oscillator(HO)-split. Let for a mo-
e ihDN needtobecomputedonlyonceandcanbereused
− mentw =1,the Hermite expansionthensuggestsa split
at each step such that the cost of the action of e ihDN
−
corresponds to N complex products. 1
A= (p2+x2), B(t)=εV (x,t)+σ(t)ψ 2, (II.22)
Fortheremainingpart,thefollowingwell-knownresult 2 I | |
is very useful
where the solution for the equation iu =Au can be ap-
′
Lemma II.1 If F is real valued, the equation proximated using a finite number of Hermite basis func-
tions, i.e.
∂
i φ(x,t)=F(x, φ(x,t))φ(x,t), (II.15) M 1
∂t | | −
φM(x,t)= cne−iEnthn(x). (II.22b)
leaves the norm invariant, φ(x,t) = φ(x,0), and then n=0
| | | | X
Since B is diagonal in the coordinate space it will act
φ(x,t)=e itF(x,φ(x,0))φ(x,0). (II.16)
− | | as a simple multiplication if we choose a number of Her-
mite basis functions and evaluate them at the points of
On the other hand, it is well known that the solutions achosenmesh(e.g. usingtheGauss-Hermitequadrature
ofthelinearSchr¨odingerequationwiththeharmonicpo- [21] or on equidistant grid points [2]). This split can be
tential of interest if all contributions from B(t) are small with
respecttoAandmethodsfornear-integrablesystemsare
i ∂ φ(x,t)= 1(p2+x2)φ(x,t) (II.17) used. Ifthefrequencyistimedependent,thecorrespond-
∂t 2 ing split is
can be expressed in terms of Hermite polynomials: 1
A(t)= (p2+ω2(t)x2), B(t)=εV (x,t)+σ(t)ψ 2.
I
2 | |
∞ (II.23)
φ(x,t)= c e iEnth (x) (II.18)
n − n In this case, it is convenient to take the time, t as a
nX=0 new coordinate as shown in (II.10). The solution of
iu =A(t)u, in the algorithm in Table I, can be approx-
where ′
imated by the Magnus expansion, e.g. (II.12) or (II.13),
En =n+ 12, hn(x)= π1/4√12nn!Hn(x)e−x2/2 bspuetctthraelsemfeatchtoodriszastiniocnesitarweonuoldt arpeqpuroirperitawtoesfoertsuosfebwaistihs
(II.19) functions and additional transformations.
and H (x) are the Hermite polynomials satisfying the
n
recursion In general, the split (i) can be considered faster and
simpler since A T can be computed in the momentum
≡
H (x)=2xH (x) 2kH , k =1,2,... space,andonecaneasilyandefficientlychangefrommo-
k+1 k k 1
− − mentum to coordinate space via FFTs. The choice (ii),
with H (x) = 1, H (x) = 2x. The weights c on the other hand, allows us to take advantage of the
0 1 n
can be computed from the initial conditions, c = structure of a near-integrable system if, roughly speak-
n
h (x)φ(x,0)dx. ing, B < A , but it requires to solve the equation
n
k k k k
R
5
for the (time-dependent) harmonic potential exactly (or In classicalmechanics, this correspondsto a kick, a drift
withhighaccuracy). Theevolutionoftheconstantoscil- and the harmonic oscillator rotation.
lator is easily computed using Hermite polynomials (see As we have seen, the exponentials can be easily com-
[2,3,21]). Asmentioned,methodstailoredforthisfamily puted by Fourier spectral methods. It is then natural to
of problems show a high performance. The main draw- ask the question if it is possible to write the solution of
back,however,isthattheevolutionfore ihaiA hastobe (II.11) as a product of exponentials which are solvable
−
carriedoutusingabasisofHermitepolynomialswhereas by spectral methods. The answer is positive and it is
theB partisadvancedinspacecoordinates[21]resulting formulated as follows.
in possibly costly basis transformations. Let us first consider the pure harmonic oscillator,
whose result was obtained in [23], and for which we
present a new proof that applies equally to the general
B. Solving the Harmonic oscillator by Fourier case.
methods
Lemma II.2 Let A = 1p2, B = 1x2 and
1 2 1 2
We propose a new method which combines the advan-
g(t)=sin(t), f(t)=tan(t/2). (II.26)
tages of both splittings. It retains the advantages of the
HO-split (ii) while being as fast to compute as the F- Then, the following property is satisfied for t <π:
split in (i). For this purpose, we briefly review some | |
basic concepts of Lie algebras. e it(A1+B1) =e if(t)A1 e ig(t)B1 e if(t)A1 (II.27)
− − − −
Given X,Y two elements of a given Lie algebra it is
=e if(t)B1 e ig(t)A1 e if(t)B1 (II.28)
well known that − − −
1 Proof. It suffices to impose that the right hand side of
eXYe−X =Y +[X,Y]+ [X,[X,Y]]+... (II.27) (or (II.28)) satisfies the equation (II.17) with the
2
identityasinitialconditionwhatcanbeverifiedbydiffer-
If X1,...,Xk arethegeneratorsofafinitedimensional entiating and simple applications of the rules (II.25a,b).
{ }
Lie algebra and A more constructive way to derive the functions f,g
makes use of the parallelism with the one-dimensional
k k
classical harmonic oscillator with Hamiltonian H =
Z = αnXn, G= anXn 1p2+ 1q2 and Hamilton equations
2 2
n=1 n=1
X X
then d q 0 1 q q
= =(A+B)
dt( p) 1 0 !( p ) ( p )
k −
G(t)=etZGe tZ = a (t)X (II.24) (II.29)
− n n
where
n=1
X
where G(t) satisfies the differential equation 0 1 0 0
A , B . (II.30)
≡ 0 0! ≡ 1 0!
G(t)=[Z,G(t)], G(0)=G. −
′
The Lie algebra generated by the matrices A,B is the
Notice that if G = x or G = p and Z = H(x,p) these
same as the Lie algebra associated to the operators
equationsareequivalenttotheequationofmotionforthe
A ,B for the Schr¨odinger equation with the harmonic
1 1
classical Hamiltonian H(x,p). It is also well-knownthat
potential (II.17).
two operators are identical on a sufficiently small time
The exact evolution operator of (II.29) is
interval if they satisfy the same first order differential
equation with the same initial conditions [22].
cos(t) sin(t)
Given X(x),P(p),W = 1(p2 + x2) and an analytic O(t)= , (II.31)
function F(x,p), we are int2erested in the following ad- −sin(t) cos(t)!
joint actions
which is an orthogonaland symplectic 2 2 matrix. For
×
e−itXF(x,p)eitX =F(x,p+tX′) (II.25a) the splitted parts, the solutions are easily computed to
e itPF(x,p)eitP =F(x tP ,p) (II.25b)
− ′ 1 f(t) 1 0
− ef(t)A = , eg(t)B =
e−itWF(x,p)eitW =F(xˆ,pˆ) (II.25c) 0 1 ! g(t) 1!
−
where X′ =dX/dx, P′ =dP/dp and and then, equating the symmetric composition
xˆ=cos(t)x+sin(t)p 1 f g 2f f2 g
efAegBefA = − · − · ,
pˆ=−sin(t)x+cos(t)p. −g 1−f ·g !
6
to (II.31), we obtain (II.26) which is valid for h π. Then, the following decomposition
| | ≤
The decomposition (II.28) is derived analogously. Using
the Baker-Campbell-Haussdorf-formula, it is clear, that e−i2t(21p2+ωL212x2)e−i2t(21p2+ωR2 21x2)
both results remain valid, up to the first singularity at
t = π, when replacing the matrices A,B by the corre- =e−if(t)12x2 e−ig(t)12p2 e−ie(t)12x2 (II.36)
±
sponding linear operators A ,B , since all computations
1 1
are done in identical Lie algebras. 2 issatisfiedfor0≤t<t∗,wheret∗ isthesmallestpositive
root of g(t).
Itisimmediatetogeneralizethisresulttotheequation The proof is similar to the previous one.
∂ 1 ω2
i ψ(x,t)= p2+µ x2 ψ(x,t),
∂t 2µ 2
(cid:18) (cid:19) III. THE HERMITE-FOURIER METHODS
µ,ω >0 by replacing (II.26) with
1 ω With the presented exact decompositions at hand, we
g = sin(ωt), f =µωtan t .
µω 2 nowsolvethediscretizedGPE(I.3)bysplittingmethods
(cid:16) (cid:17) using the symmetric compositions (II.3) and (II.4). Let
This result is valid for t <t π/ω.
| | ∗ ≡ us first consider the case w = 1 and take the HO split
The following theorems extend this idea to decompo-
A=A +B
sitions of operators that appear after the approximation 1 1
of the time dependent parts via (II.12) or by the compo-
sition (II.13). ψh[2,]A =e−ih(A1+B1)/2e−ihBe−ih(A1+B1)/2 (III.1)
ψ[2] =e ihB/2e ih(A1+B1)e ihB/2. (III.2)
Theorem II.3 Let α,β,γ be constants, η = αγ β2 h,B − − −
−
and
p Replacing the exponentials e ih(A1+B1) by (II.27) or
−
g(t)=γ/η sin(ηt), (II.32)
(II.28),weobtainfourdifferentmethodswhosecomputa-
·
1 β tional costs differ considerably. At first sight, using the
f(t)= 1 cos(ηt)+ sin(ηt) ,
g(t) − η FSAL property, both (III.1) and (III.2) are equivalent
(cid:18) (cid:19)
from the computational point of view and require one
1 β
e(t)= 1 cos(ηt) sin(ηt) . exponential of B and another one of A +B per step.
1 1
g(t) − − η
(cid:18) (cid:19) However, a significant difference arises when we plug in
Then, the following decomposition holds for 0 t<π/η: the decompositions (II.27) or (II.28). Only the combi-
≤
nation (III.2) with (II.28) yields a method that involves
e−i2t(αx2+β(xp+px)+γp2) only a single FFT call per step
=e−if(t)12x2 e−ig(t)21p2 e−ie(t)12x2. (II.33) ψ[2] =e ihB/2e ih(A1+B1)e ihB/2
h − − −
Proof. The proof follows the lines of the proof of
=e ihB/2e if(h)B1e ig(h)A1e if(h)B1e ihB/2
Lemma II.2. The evolution operator associated to the − − − − −
classicalHamiltonianH = 1 αx2+2βxp+γp2 isgiven =e i(hB/2+f(h)B1)e ig(h)A1e i(hB/2+f(h)B1)
2 − − −
by (III.3)
(cid:0) (cid:1)
cos(ηt)+ β sin(ηt) γ sin(ηt)
−αη sinη(ηt) cos(ηtη)− βη sin(ηt)!. (II.34) aFnodr asnolyveostheexracctolmybthineathiaornm, omnoicreokscinilelatticorteformr s|hh|a<veht∗o.
becomputedperstepsincetheFSALpropertycannotbe
andequalitytotherighthandsideof (II.33)isverifiedby
exploited to full extent and hence result in more costly
straight-forwardcomputationofthematrixexponentials.
[24] methods for the same accuracy.
Thesolutionisvaliduntilthefirstsingularityatt=π/η.
The generalcomposition(II.6) canbe rewrittenin the
Using (II.25a,b), it can be checked that both sides of
same way by replacing each flow e iaiA in (II.6) by the
(II.33) satisfy the same differential equation and initial −
conditions. 2 composition (II.28)
Tsinh(eωorte/m2) fIoIr.4kL=etLω,Rk a∈ndR, ck = cos(ωkt/2), sk = Φh ≡e−i(hbmB+αmB1)e−ig(amh)A1e−i(hbmB+αm−1B1)
k e i(hb1B+α1B1)e ig(a1h)A1e iα0B1 (III.4)
− − −
g(t)=s c /ω +c s /ω , (II.35) ···
L R L L R R
1 ωL where αk =f(ak+1h)+f(akh), k=0,1,...,m+1 with
f(t)= g(t) 1−cLcR+ ω sLsR , a0 =am+1 =0. Thismethodisvalidfor aih <h∗, i=
(cid:18) R (cid:19) 1,...,m and requires only m calls of th|e F|FT and its
1 ω
e(t)= 1 c c + Rs s . inverse,but reaches the same accuracy as if the Hermite
L R L R
g(t)(cid:18) − ωL (cid:19) functions were used.
7
In the more general case with a time-dependent fre-
0 0
quency, ω(t), starting from the HO splitting, the time-
dependent part is first approximated by Magnus expan-
−5 −5
sions (II.12) or (II.13). Theorems II.3 and II.4 then pro-
videdecompositionstowritetheproductofexponentials −10 −10
in a similar way as in (III.4) but now α = f(a h)+
k k+1
e(akh) and it is valid for aih <t∗, i=1,...,m where −1−5π −2 0 2 π −153 3.05 3.1 π
| |
t is the first zero of g(t). At each stage, one has to
∗
compute f(a h),g(a h),e(a h) from ω(t). FIG. 1. Error in logarithmic scale for the integration of the
i i i
As in the previous case, it only requires m calls of groundstateoftheHarmonicpotentialusingthesplit(II.28)
for T ∈ [−π,π] (integration forward and backward in time).
the FFT and its inverse, like the standard Fourier pseu-
The left and right panels show the 2-norm error and a zoom
dospectral methods. For stability reasons, it seems
about T =π, respectively.
convenient to look for splitting method whose value of
max a is as small as possible.
i i
{| |}
1,...,M, j =1,...,N =512withh (x)givenin(II.18),
n
IV. NUMERICAL EXAMPLES x [ 10,10]. For δ = 1 , round off accuracyis achieved
∈ − 10
with M = 8 while for δ = 2 it is necessary to take
We analyze first the performance of the methods con- M = 29. We observe that the Hermite decomposition
sideredinthisworkfortheone-dimensionalproblem(I.1) is very sensitive to the initial conditions [21]. The Her-
with µ = ω2 = 1, and the pure harmonic trap, i.e. mitebasisworksefficientlyasfarastheinitialconditions
V =0. aswellastheevolutionthereofcanbeaccuratelyapprox-
I
To illustrate the validity of the decomposition pre- imatedusingafewnumberofbasiselementsandonehas
sented in Lemma II.2, we first consider the linear prob- tokeepinmindthat,fornonlinearproblems,thenumber
lem (σ =0). We take the groundstate att=0 as initial ofbasisfunctionsnecessarytoreachagivenaccuracycan
condition whose exact solution is given by vary along the time integration.
Next, we study the following values for the nonlinear-
ψ(x,t)= 1 e−it/2e−x2/2. ity parameter: σ = 10−2,1,102. The case σ = 10−2
π1/4√2nn! illustrates the performance of the new methods if ap-
plied to problems (linear or nonlinear) which are small
Wediscretizeontheinterval[ 10,10],toensurethewave
− perturbationsoftheHarmonicpotentialwhereastheval-
function and its first derivatives vanish up to round off
ues σ =1,102 are large enough to demonstrate the non-
at the boundaries, and sample it at N = 1024 equidis-
linearity effects on the approximation properties of the
tant grid points. We integrate, with only one time step,
methods. Physically, σ is proportional to the number of
from t = 0 to T for T [ π,π], i.e. forward and back-
∈ − particles in a Bose-Einsteincondensate and to the inter-
ward in time. We measure the integrated error in the
action strength [5].
wave function, u (T) u (T) , where u (T) denotes
k ex − a k2 a For all cases, we choose the initial condition ψ(x,0)=
the approximate numerical solution obtained using the
ρe (x 1)2/2, with ρ a normalizing constant. We show in
split (II.28) and u (T) is the exact solution at the dis- − −
ex Fig. 2 the value of ψ(x,t)2 at the initial and final time.
cretizedmesh. Theresultofthiscomparisonisillustrated | |
The spatial interval is adjusted to ensure the wave func-
in Fig. 1 (left). The split (II.28) reproduces,for T <π,
| | tion vanishes (up to round off) at the boundaries, here
the exact solution up to round off, as expected. The
[ 20,20] for σ = 0.01 and [ 30,30] for both σ = 1 and
right panel in Fig. 1 displays a zoom near a singularity
− −
σ = 100 where the wave function moves faster (we only
where the error grows rapidly due to double precision
show the interval x [ 5,5]). One can appreciate that
arithmetic.
∈ −
for strong nonlinearities the wave function can consid-
We analyze how the approximation properties of the
erably penetrate the potential barrier and one expects
Hermite decomposition (II.22b) strongly depend on the
that an accurate approximation of these wave functions
function in question and on the chosen number of ba-
requiresalargenumber ofHermite functions whenusing
sis functions, M. We compute the M required to
(IV.1) and hence renders this procedure inappropriate.
reach round-off precision for the evolution of a dis-
placed ground state as initial condition, ψ (x,0) =
δ
e (x δ)2/2/(π1/4√2nn!) from t = 0 to T = 10 in one
− − A. HO-split versus F-split
time step. From initial conditions computed on a mesh,
this can be accomplished as follows [2]:
We analyzenowthe advantagesofthe HO-splitversus
u (T)=e iT(A1+B1)u KTe iTD1Ku (IV.1) the F-split as given in (II.23) and (II.21).
ex − 0 − 0
∼ In a first experiment, we fix the symmetric second or-
where D = diag 1,3,...,2M 1 , M is the number of derBAB composition(III.2)andapplyitforbothsplits.
1 {2 2 2− }
basis elements considered, and K = h (x ), i = For the HO-split, we compute the harmonic part either
i,j i 1 j
−
8
σ=0.01, T=100 σ=1, T=100 σ=100, T=10
TABLE II.Coefficients for several splitting methods.
0.6 0.6 0.6
The 6-stage 4th-order: SRKN 4
0.4 0.4 0.4 6
2Ψ| b =0.0829844064174052 a =0.245298957184271
| 0.2 0.2 0.2 1 1
b =0.396309801498368 a =0.604872665711080
2 2
0 0 0 b =−0.0390563049223486 a =1/2−(a +a )
−5 0 5 −5 0 5 −5 0 5 3 3 1 2
x x x b =1−2(b +b +b ) a =a , a =a , a =a
(N=512) (N=1024) (N=1024) 4 1 2 3 4 3 5 2 6 1
b =b , b =b , b =b
5 3 6 2 7 1
The4-stage (8,2) method: NI
FIG. 2. Exact evolution at t=T (solid line) from the initial 4
conditions given by ψ(x,0)=ρe−(x−1)2/2 (dashed line). The b1 =1/20 a1=1/2−p3/28
numberof grid points is given byN. b =49/18 a =1/2−a
2 2 1
b =1−2(b +b ) a =a , a =a
3 1 2 3 2 4 1
b =b , b =b
4 2 5 1
with the decomposition (II.28) or in the Hermite basis The5-stage (8,4) method: NI5
(IV.1) with different numbers of basis terms. b =0.811862738544516 a =−0.00758691311877447
1 1
The parameters, initial conditions and final times are b =−0.677480399532169 a =0.317218277973169
2 2
taken as previously for Fig. 2. Given a splitting method
b =1/2−(b +b ) a =1−2(a +a )
3 1 2 3 1 2
X,wedenotebyX ,X andX M itsimplementations
F H H b =b , b =b , b =b a =a , a =a
with the F-split, the HO-split using the Hermite-Fourier 4 3 5 2 6 1 4 2 5 1
method and the HO-split using M Hermite basis func-
tions in (IV.1), respectively. We measure the error ver-
sus the number of exponentials which can be considered ods, LF, are compared with the second order NI4(8,2)
proportional to the computational cost and plot the re- methods. In the lower panels we compare the RKN64
sults in Fig. 3. As expected, HO-splits are advantagous methods against the (8,4) methods jointly with the best
if the system is close to a harmonic oscillator, i.e. for among the previous second order methods.
σ = 0.01,1, and if the initial conditions are accurately For a weak nonlinearity, when the system can be con-
approximated by a few terms in the Hermite expansion. sideredasaperturbedharmonicoscillator,weclearlyob-
Ontheotherhand,forstrongnonlinearitiesσ =100,the servethat the HO-split is superior to the F-split. In this
Hermite polynomial based HO split shows a very lim- case,witharelativelysmallnumberofHermitefunctions,
ited performance, c.f. the large number of basis terms it is possible to approach accurately the solution, but
in the right panel of Fig. 3. We stress that if this tech- this procedure has a limited accuracywhich candeterio-
niqueisusedfornonlinearitiesthenumberofbasisterms ratealongthetimeintegrationanddependsontheinitial
should be increased along the time integration and fix- conditions. In addition, the methods addressed to per-
ing it bounds the maximally achievable accuracy and its turbedproblemsshowthebestperformance: The(8,2)
H
limitdependsontheinitialconditionandthestrengthof method performs best among the compared when a rel-
the nonlinearity. atively low accuracy is desired and the (8,4) method
H
The Hermite-Fourier method proposed in this work takes its place for higher accuracies.
(using the composition (II.28)) is clearly superior for Figure 5 shows the results for σ = 1. It is qualita-
weak perturbations and it keeps similar performance to tively similar to the previous case yet the HO-split does
the F-split for strong nonlinearities. not outperform the plain F-split (II.21) as significantly
Finally,weanalyzetheperformanceofdifferenthigher as before. Nevertheless, it is important to observe that,
order splitting methods which are useful when high ac- again,thebestresultisobtainedfortheHO-split. Notice
curacies are desired. The following methods (whose co- thatahighernumberofHermitebasisfunctionsisneces-
efficients are collected in Table II for the convenience of sarytoachievethesameaccuracyastheHermite-Fourier
the reader) are considered: decomposition.
– RKN64 (the 6-stage fourth-order method from [9]). Figure 6 shows the results for σ = 100. The HO-split
Thisis aPartitionedRunge-Kutta-Nystr¨ommethodand cannotbe expected to be particularlyuseful because the
it is designed for the case where [B,[B,[B,A]]] = 0, be- system is far from being a harmonic oscillator. From
ing the case for both the F-split and the HO-split. Fig. 2, we expect a great number of Hermite basis func-
– NI4(8,2) (the 4-stage (8,2) BAB method from [11]). tionstobe requiredforasufficientlyaccurateexpansion.
Thismethodisaddressedtoperturbedsystems. Oneex- The results in Fig. 6 demonstrate this rather intuitive
pects a high performance if the contribution from B is expectation, i.e. almost negligible precision despite the
small. large number of basis terms M = 150. Remarkably, the
– NI5(8,4) (the 5-stage (8,4) BAB method from [11]). proposedHO decomposition does not show these limita-
We analyze in Figures 4-6 the three problems speci- tions and reaches the precision of the F-split (II.21) be-
fied in Fig. 3. In the upper panels, the Leap-Frog meth- causewearesolvingtheharmonicpotentialexactlyupto
9
σ=0.01, T=100, N=512 σ=1, T=100, N=1024 σ=100, T=10, N=1024
1 0.5 0.5
0 0 0
−0.5
−1 −0.5
T)||2 T)||2 T)||2 −1
φ−( −2 φ−( −1 φ−( −1.5
T) T) T)
Ψ( Ψ( Ψ(
g||10−3 g||10−1.5 g||10 −2
lo lo lo −2.5
−4 LF −2 LF LF
F F F
LF LF −3 LF
H H H
−5 LF 10 −2.5 LF 10 LF 100
H H −3.5 H
LF 30 LF 30 LF 150
H H H
−6 −3 −4
2 2.5 3 3.5 4 2.5 3 3.5 4 3 3.2 3.4 3.6 3.8 4
log (Nr of Exp) log (Nr of Exp) log (Nr of Exp)
10 10 10
FIG. 3. Error versusthenumberof exponentials in logarithmic scale for different splittings for the Leap-Frog method.
spectral accuracy. For this problem, we observethat the equivalent to take in (IV.4) t = a h and the parameters
i
(8,2)methodhasthebestperformancewhenarelatively α,β,γ are given by:
lowaccuracyisdesired,the(8,4)methodshowsthe best
α= 1 (5ω +8ω +5ω )
ptheerfomremthaondceoffocrhmoiceedifuomr haigchcuerraacciecsuraancdiest.he RKN64 is +(1a8ih)2(117(ω22+ω2)+3 8ω2+ω ω +ω ω 37ω ω ),
486 4 1 3 2 1 2 2 3− 2 1 3
β =a h 5(ω ω )( 1 + (aih)2(5ω +8ω +5ω )),
i 3 3− 1 12 3240 1 2 3
B. Time-dependent harmonic oscillator perturbed γ =1+q(aih)2(ω 2ω +ω )
54 1− 2 3
by weak quartic anharmonicity
withω =ω(t +c h), i=1,2,3andc =1/2 √15/10,
i n i 1
−
c = 1/2, c = 1/2+√15/10, corresponding to a sixth-
We consider now a harmonic oscillator with time- 2 3
order Gaussian quadrature rule. The obtained operator
dependent frequency and perturbed by a weak static
is then decomposed according to Theorem II.3.
quartic anharmonicity
The results are given in Fig. 7 and corroborate the
∂ 1 1 1 superiority of the HO split.
i ψ = p2+ ω2(t)x2 ψ+ε x4ψ. (IV.2)
∂t 2 2 Q4 Another interesting example is given by an intense
(cid:18) (cid:19) short pulse modeled via
We first consider the case ω2(t) = A(1+ǫcos(wt)) with
At
w = 1/2, A = 4, ǫ = 0.1, ε = 0.01. As reference, we ω(t)=w 1+ .
Q 0 cosh2(B(t 2))
take a highly accurate numerical approximation as ex- (cid:18) − (cid:19)
act solution and restrict the spatial domain to [ 20,20] Varying the parameters A and B, the pulse can be
−
for all experiments in this subsection . We compare sharpened while keeping its time-average and hence its
the Hermite-Fouriermethodwith the plainFouriersplit, strength relative to the anharmonicity constant. Fig-
sinceHermitepolynomialsarenotappropriateinatime- ure 8 shows the results obtained for a relatively slow
dependentsetting. Forfastoscillatingsystemsandifhigh variation of the harmonic potential, for the parameters:
accuracyis needed,the two-exponentialfourth-orderap- w = 4, A =0.25, B = 2. Again, the advantageousness
0
proximationoftheharmonicoscillator(II.13)canbeim- of the presented decomposition can be appreciated. It is
proved by taking, for example, a higher order Magnus alreadynoticeable,thattheerrorintroducedbythetime-
expansion (II.12). As we have seen, the solution of dependence becomes dominant, and this effect increases
for more rapidly varying potentials, e.g. for B 1. In
1 1 ≫
iU′ = p2+ ω2(t)x2 U (IV.3) that case, higher order approximations of the Magnus
2 2 expansion are necessary to maintain the benefits of the
(cid:18) (cid:19)
Hermite composition.
can be written as
U(t,0)=e−i2t(αx2+β(xp+px)+γp2), (IV.4)
V. CONCLUSIONS
andwe have considered,for example, a sixth-orderMag-
nusintegrator[17]toapproximatetheevolutionoperator We have presented new Fourier methods for the nu-
foronefractionaltimestep,a h,i.e. U(t,t+a h). Thisis merical integration of perturbations to the time depen-
i i
10
monicpotential. Themethodspresentedinthisworkex-
σ=0.01, T=100, N=512
0 σ=1, T=100, N=1024
LFF
LF 0
−1 H
LF 30
H −0.5
−2 (8,2)
F
Ψφlog||(T)−(T)||102 −−−543 ((88,,22))HH 30 Ψφ||(T)−(T)||102 −−21−−..2155 LFF
g
−6 lo −3 LFH
LF 30
H
−7 −3.5 (8,2)
F
(8,2)
−8 −4 H
(8,2) 30
2.5 3 3.5 4 4.5 H
log (Nr of Exp) −4.5
10 2.5 3 3.5 4 4.5
σ=0.01, T=100, N=512 log (Nr of Exp)
0 σ=1, T10=100, N=1024

−1 0
−2 −1
−3 −2
Ψφ||(T)−(T)||102 −−−654 ((88,,24))HF Ψφ(T)−(T)||2 −−43 (8,2)H
log −7 (8,4)H ||10 −5 (8,4)F
g
(8,4) 50 o (8,4)
−8 H l −6 H
RKN4 (8,4) 80
6 F H
−9 RKN4 −7 RKN4
6 H 6 F
−10 RKN4 50 RKN4
6 H −8 6 H
RKN4 80
2.5 3 3.5 4 6 H
log (Nr of Exp) −9
10 2.5 3 3.5 4 4.5
log (Nr of Exp)
10
FIG.4. Comparisonofsecondorder(upperpanel)andfourth
order (lower panel) methods for the different splittings and FIG. 5. Same as Fig. 4 for σ=1.
decompositions discussed in thetext and σ=0.01.
tendtoperturbedharmonicpotentialsinlinearquantum
dent harmonic oscillator which are useful for both both
mechanics, c.f. section IVB, where it is straightforward
the Gross-Pitaevskii equation as well as for the linear
to generalize the results to higher dimensions.
Schr¨odinger equation. Fourier methods have shown a
high performance in solving many different problems
which can be split into the kinetic part and a remain-
ACKNOWLEDGMENTS
der that is diagonal in the coordinate space. We have
extended the Fourier methods to perturbations of the
time-dependent harmonic potential, and refer to them We would like to acknowledge the referees for their
as Hermite-Fourier methods. They solve the linear careful reading and suggestions. The authors acknowl-
Schr¨odinger equation with a time-dependent harmonic edge the support of the Generalitat Valenciana through
potential to the desired order using corresponding Mag- theprojectGV/2009/032. TheworkofS.Blaneshasalso
nus expansions andup to the accuracygivenby the spa- beenpartiallysupportedbyMinisteriodeCienciaeInno-
tial discretization. These methods are fast to compute vacio´n(Spain)underthecoordinatedprojectMTM2010-
since FFTs can be applied and show a high accuracy 18246-C03 (co-financed by the ERDF of the European
when the problem is a small perturbation of the har- Union).