quspin.tools.evolution.evolve

quspin.tools.evolution.evolve(v0, t0, times, f, solver_name='dop853', real=False, stack_state=False, verbose=False, imag_time=False, iterate=False, f_params=(), **solver_args)[source]

Implements (imaginary) time evolution for a user-defined first-order ODE.

The function can be used to study nonlinear semiclassical dynamics. It can also serve as a pre-configured ODE solver in python, without any relation to other QuSpin objects.

Parameters:
v0numpy.ndarray

Initial state.

t0float

Initial time.

timesnumpy.ndarray

Vector of times to compute the time-evolved state at.

ffunction

User-defined function to solve first-order ODE (see Examples):

\[v'(t) = f(v(t),t)\qquad v(t_0)=v_0\]
f_paramstuple, optional

A tuple to pass all parameters of the function f to ODE solver. Default is f_params=().

iteratebool, optional

If set to True, creates a generator object for the time-evolved the state. Default is False.

solver_namestr, optional

Scipy solver integrator name. Default is dop853.

See scipy integrator (solver) for other options.

solver_argsdict, optional

Dictionary with additional scipy integrator (solver) arguments.

realbool, optional

Flag to determine if f is real or complex-valued. Default is False.

imag_timebool, optional

Must be set to True when f defines imaginary-time evolution, in order to normalise the state at each time in times. Default is False.

stack_statebool, optional

If f is written to take care of real and imaginary parts separately (see Examples), this flag will return a single complex-valued solution instead of the real and imaginary parts separately. Default is False.

verbosebool, optional

If set to True, prints normalisation of state at teach time in times.

Returns:
obj
Can be either one of the following:

  • numpy.ndarray containing evolved state against time.

  • generator object for time-evolved state (requires iterate = True).

Examples

The following example shows how to use the evolve() function to solve the periodically-driven Gross-Pitaevskii equation (GPE) on a one-imensional lattice. The GPE has a linear part, comprising the kinetic energy and the external potentials (e.g. a harmonic trap), and a nonlinear part which describes the interactions.

Below, in a few steps we show how to use the functionality of the evolve() function to solve the GPE on a one-dimensional lattice for periodically-driven particles in a harmonic trap:

\[i\dot\varphi_j(t) = -J\left( e^{-iA\sin\Omega t}\varphi_{j-1}(t) + e^{+iA\sin\Omega t}\varphi_{j+1}(t) \right) + \kappa_\mathrm{trap}\varphi_j(t) + U|\varphi_j(t)|^2\varphi_j(t)\]

where \(j\) labels the lattice sites, \(J\) is the lattice hopping amplitude, \(\kappa_\mathrm{trap}\) is the strength of the harmonic trap, and \(U\) – the mean-field interaction.

Let us start by defining the single-particle Hamiltonian \(H(t)\).

 1from quspin.basis import boson_basis_1d  # Hilbert space spin basis
 2from quspin.tools.evolution import evolve  # ODE evolve tool
 3from quspin.tools.Floquet import Floquet_t_vec  # stroboscopic time vector
 4import numpy as np  # generic math functions
 5from six import iteritems  # loop over elements of dictionary
 6
 7#
 8L = 50  # number of lattice sites
 9i_CM = L // 2 - 0.5  # centre of chain
10#
11### static model parameters
12J = 1.0  # hopping
13kappa_trap = 0.002  # harmonic trap strength
14U = 1.0  # mean-field (GPE) interaction
15#
16### periodic driving
17A = 1.0  # drive amplitude
18Omega = 10.0  # drive frequency
19
20
21def drive(t, Omega):
22    return np.exp(-1j * A * np.sin(Omega * t))
23
24
25def drive_conj(t, Omega):
26    return np.exp(+1j * A * np.sin(Omega * t))
27
28
29drive_args = [Omega]  # drive arguments
30t = Floquet_t_vec(Omega, 30, len_T=1)  # time vector, 30 stroboscopic periods
31#
32### site-couping lists
33hopping = [[-J, i, (i + 1) % L] for i in range(L)]
34trap = [[kappa_trap * (i - i_CM) ** 2, i] for i in range(L)]
35#
36### operator strings for single-particle Hamiltonian
37static = [["n", trap]]
38dynamic = [["+-", hopping, drive, drive_args], ["-+", hopping, drive_conj, drive_args]]
39# define single-particle basis
40basis = boson_basis_1d(
41    L, Nb=1, sps=2

Next, we define the GPE

\[-i\dot\varphi(t) = H(t)\varphi(t) + U |\varphi(t)|^2 \varphi(t)\]

and solve it using evolve():

43#
44### build Hamiltonian
45H = hamiltonian(static, dynamic, basis=basis, dtype=np.complex128)
46# calculate eigenvalues and eigenvectors of free particle in a harmonic trap
47E, V = H.eigh(time=0)
48# initial state normalised to one partcile per site
49phi0 = V[:, 0] * np.sqrt(L)
50
51
52#######
53def GPE(time, phi):
54    """Solves the complex-valued time-dependent Gross-Pitaevskii equation:"""
55    # integrate static part of GPE

Since the GPE is a complex-valued equation, the above code requires the use of a complex-valued ODE solver [which is done by evolve() under the hood, so long as no solver is explicitly specified].

An alternative way to solve the GPE using a real-valued solver might be useful to speed-up the computation. This can be achieved by decomposing the condensate wave function into a real and imaginary part, and proceeds as follows:

The goal is to solve:

\[-i\dot\varphi(t) = H(t)\varphi(t) + U |\varphi(t)|^2 \varphi(t)\]

for the complex-valued \(\varphi(t)\) by re-writing it as a real-valued vector phi=[u,v] where \(\varphi(t) = u(t) + iv(t)\). The real and imaginary parts, \(u(t)\) and \(v(t)\), have the same array dimension as \(\phi(t)\).

In the most general form, the single-particle Hamiltonian can be decomposed as \(H(t)= H_{stat} + f(t)H_{dyn}\), with a complex-valued driving function \(f(t)\). Then, the GPE can be cast in the following real-valued form:

\[\dot u(t) = +\left[H_{stat} + U\left(|u(t)|^2 + |v(t)|^2\right) \right]v(t) + Re[f(t)]H_{dyn}v(t) + Im[f(t)]H_{dyn}u(t)\]
\[\dot v(t) = -\left[H_{stat} + U\left(|u(t)|^2 + |v(t)|^2\right) \right]u(t) - Re[f(t)]H_{dyn}u(t) + Im[f(t)]H_{dyn}v(t)\]
 59    # integrate dynamic part of GPE
 60    for fun, Hdyn in iteritems(H.dynamic):
 61        phi_dot += -1j * fun(time) * Hdyn.dot(phi)
 62
 63    return phi_dot
 64
 65
 66# solving the complex-valued GPE takes one line
 67phi_t = evolve(phi0, t.i, t.vals, GPE)
 68
 69
 70########
 71def GPE_real(time, phi, H, U):
 72    """Solves the Gross-Pitaevskii equation, cast into real-valued form so it can be solved with a
 73    real-valued ODE solver.
 74
 75    """
 76    # preallocate memory for phi_dot
 77    phi_dot = np.zeros_like(phi)
 78    # read off number of lattice sites (array dimension of phi)
 79    Ns = H.Ns
 80    # static single-particle part
 81    phi_dot[:Ns] = H.static.dot(phi[Ns:]).real
 82    phi_dot[Ns:] = -H.static.dot(phi[:Ns]).real
 83    # static GPE interaction
 84    phi_dot_2 = np.abs(phi[:Ns]) ** 2 + np.abs(phi[Ns:]) ** 2
 85    phi_dot[:Ns] += U * phi_dot_2 * phi[Ns:]
 86    phi_dot[Ns:] -= U * phi_dot_2 * phi[:Ns]
 87    # dynamic single-particle term
 88    for func, Hdyn in iteritems(H.dynamic):
 89        fun = func(time)  # evaluate drive
 90        phi_dot[:Ns] += (
 91            +(fun.real) * Hdyn.dot(phi[Ns:]) + (fun.imag) * Hdyn.dot(phi[:Ns])
 92        ).real
 93        phi_dot[Ns:] += (
 94            -(fun.real) * Hdyn.dot(phi[:Ns]) + (fun.imag) * Hdyn.dot(phi[Ns:])
 95        ).real
 96
 97    return phi_dot
 98
 99
100# define ODE solver parameters
101GPE_params = (H, U)
102# solving the real-valued GPE takes one line
103phi_t = evolve(phi0, t.i, t.vals, GPE_real, stack_state=True, f_params=GPE_params)

The flag stack_state=True is required for evolve() to handle the complex-valued initial condition properly, as well as to put together the output solution as a complex-valued vector in the end. Since the real-valued ODE solver allows to parse ODE parameters, we can include them in the user-defined ODE function and use the flag f_params. Notice the elegant way python allows one to circumvent the usage of this variable in the complex-valued example above.