Exact Solutions

Exact Solutions#

Here we present the various exact solutions described in gpe/exact_solutions.py

Harmonic Oscillator#

Here we construct an exact solution for a trapped BEC in a harmonic oscillator:

\[\begin{gather*} \psi_0(x) = \sqrt{n_0}e^{-x^2/2\sigma^2}, \qquad V(x) = \frac{m\omega^2x^2}{2} - \frac{\hbar^2}{2m\sigma^2} - gn_0e^{-x^2/\sigma^2}, \qquad \sigma = \sqrt{\frac{\hbar}{m\omega}}. \end{gather*}\]

This potential has the non-linear piece subtracted so that the constructed solution is an exact solution with zero chemical potential (i.e. a stationary state) and energy:

\[\begin{gather*} E = \int_{-\infty}^{\infty}\d{x}\;\frac{-gn_0^2}{2}e^{-2x^2/\sigma^2} = -\sqrt{\frac{\pi}{2}}\frac{gn_0^2\sigma}{2} \end{gather*}\]

%pylab inline --no-import-all
from gpe.imports import *
from gpe.minimize import MinimizeState
from gpe.exact_solutions import HarmonicOscillator

%pylab is deprecated, use %matplotlib inline and import the required libraries.
Populating the interactive namespace from numpy and matplotlib
[I 20:00:06 numexpr.utils] NumExpr defaulting to 2 threads.

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[2], line 2
get_ipython().run_line_magic('pylab', 'inline --no-import-all')
----> 2 from gpe.imports import *
from gpe.minimize import MinimizeState
from gpe.exact_solutions import HarmonicOscillator

File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/imports.py:28
from mmfutils.plot import imcontourf  # noqa: E402
from gpe.minimize import MinimizeState  # noqa: E402
---> 28 from gpe.utils import evolve_to, evolve, evolves  # noqa: E402
from gpe.plot_utils import MPLGrid  # noqa: E402
from pytimeode.evolvers import EvolverSplit, EvolverABM  # noqa: E402

File ~/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/utils.py:35
from persist.objects import Archivable
from persist.archive import Archive
---> 35 from pytimeode.evolvers import EvolverABM
from pytimeode.mixins import ArrayStateMixin
from pytimeode.interfaces import implementer, IStateForABMEvolvers

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/pytimeode/__init__.py:2
from . import interfaces
----> 2 from . import mixins
from . import evolvers
__all__ = ["interfaces", "mixins", "evolvers"]

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/pytimeode/mixins.py:39
from .interfaces import (
   IState,
   IStateApply,
   (...)     36     implementer,
)
from . import interfaces
---> 39 from .utils import expr
__all__ = [
   "StateMixin",
   "StatesMixin",
   (...)     47     "ArraysStateWithBraketMixin",
]
@implementer(IState)
class StateMixin:

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/pytimeode/utils/expr.py:18
import numpy as np
try:
---> 18     import sympy
except ImportError:
   sympy = None

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/__init__.py:77
from .logic import (to_cnf, to_dnf, to_nnf, And, Or, Not, Xor, Nand, Nor,
       Implies, Equivalent, ITE, POSform, SOPform, simplify_logic, bool_map,
       true, false, satisfiable)
from .assumptions import (AppliedPredicate, Predicate, AssumptionsContext,
       assuming, Q, ask, register_handler, remove_handler, refine)
---> 77 from .polys import (Poly, PurePoly, poly_from_expr, parallel_poly_from_expr,
       degree, total_degree, degree_list, LC, LM, LT, pdiv, prem, pquo,
       pexquo, div, rem, quo, exquo, half_gcdex, gcdex, invert,
       subresultants, resultant, discriminant, cofactors, gcd_list, gcd,
       lcm_list, lcm, terms_gcd, trunc, monic, content, primitive, compose,
       decompose, sturm, gff_list, gff, sqf_norm, sqf_part, sqf_list, sqf,
       factor_list, factor, intervals, refine_root, count_roots, all_roots,
       real_roots, nroots, ground_roots, nth_power_roots_poly, cancel,
       reduced, groebner, is_zero_dimensional, GroebnerBasis, poly,
       symmetrize, horner, interpolate, rational_interpolate, viete, together,
       BasePolynomialError, ExactQuotientFailed, PolynomialDivisionFailed,
       OperationNotSupported, HeuristicGCDFailed, HomomorphismFailed,
       IsomorphismFailed, ExtraneousFactors, EvaluationFailed,
       RefinementFailed, CoercionFailed, NotInvertible, NotReversible,
       NotAlgebraic, DomainError, PolynomialError, UnificationFailed,
       GeneratorsError, GeneratorsNeeded, ComputationFailed,
       UnivariatePolynomialError, MultivariatePolynomialError,
       PolificationFailed, OptionError, FlagError, minpoly,
       minimal_polynomial, primitive_element, field_isomorphism,
       to_number_field, isolate, round_two, prime_decomp, prime_valuation,
       galois_group, itermonomials, Monomial, lex, grlex,
       grevlex, ilex, igrlex, igrevlex, CRootOf, rootof, RootOf,
       ComplexRootOf, RootSum, roots, Domain, FiniteField, IntegerRing,
       RationalField, RealField, ComplexField, PythonFiniteField,
       GMPYFiniteField, PythonIntegerRing, GMPYIntegerRing, PythonRational,
       GMPYRationalField, AlgebraicField, PolynomialRing, FractionField,
       ExpressionDomain, FF_python, FF_gmpy, ZZ_python, ZZ_gmpy, QQ_python,
       QQ_gmpy, GF, FF, ZZ, QQ, ZZ_I, QQ_I, RR, CC, EX, EXRAW,
       construct_domain, swinnerton_dyer_poly, cyclotomic_poly,
       symmetric_poly, random_poly, interpolating_poly, jacobi_poly,
       chebyshevt_poly, chebyshevu_poly, hermite_poly, hermite_prob_poly,
       legendre_poly, laguerre_poly, apart, apart_list, assemble_partfrac_list,
       Options, ring, xring, vring, sring, field, xfield, vfield, sfield)
from .series import (Order, O, limit, Limit, gruntz, series, approximants,
       residue, EmptySequence, SeqPer, SeqFormula, sequence, SeqAdd, SeqMul,
       fourier_series, fps, difference_delta, limit_seq)
from .functions import (factorial, factorial2, rf, ff, binomial,
       RisingFactorial, FallingFactorial, subfactorial, carmichael,
       fibonacci, lucas, motzkin, tribonacci, harmonic, bernoulli, bell, euler,
   (...)    138         Znm, elliptic_k, elliptic_f, elliptic_e, elliptic_pi, beta, mathieus,
       mathieuc, mathieusprime, mathieucprime, riemann_xi, betainc, betainc_regularized)

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/polys/__init__.py:79
__all__ = [
   'Poly', 'PurePoly', 'poly_from_expr', 'parallel_poly_from_expr', 'degree',
   'total_degree', 'degree_list', 'LC', 'LM', 'LT', 'pdiv', 'prem', 'pquo',
   (...)     65     'field', 'xfield', 'vfield', 'sfield'
]
from .polytools import (Poly, PurePoly, poly_from_expr,
       parallel_poly_from_expr, degree, total_degree, degree_list, LC, LM,
       LT, pdiv, prem, pquo, pexquo, div, rem, quo, exquo, half_gcdex, gcdex,
   (...)     76         nth_power_roots_poly, cancel, reduced, groebner, is_zero_dimensional,
       GroebnerBasis, poly)
---> 79 from .polyfuncs import (symmetrize, horner, interpolate,
       rational_interpolate, viete)
from .rationaltools import together
from .polyerrors import (BasePolynomialError, ExactQuotientFailed,
       PolynomialDivisionFailed, OperationNotSupported, HeuristicGCDFailed,
       HomomorphismFailed, IsomorphismFailed, ExtraneousFactors,
   (...)     91         MultivariatePolynomialError, PolificationFailed, OptionError,
       FlagError)

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/polys/polyfuncs.py:10
from sympy.polys.polyoptions import allowed_flags, build_options
from sympy.polys.polytools import poly_from_expr, Poly
---> 10 from sympy.polys.specialpolys import (
   symmetric_poly, interpolating_poly)
from sympy.polys.rings import sring
from sympy.utilities import numbered_symbols, take, public

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/polys/specialpolys.py:298
   return dmp_mul(f, h, n, K), dmp_mul(g, h, n, K), h
# A few useful polynomials from Wang's paper ('78).
--> 298 from sympy.polys.rings import ring
def _f_0():
   R, x, y, z = ring("x,y,z", ZZ)

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/polys/rings.py:31
from sympy.polys.polyoptions import (Domain as DomainOpt,
                                    Order as OrderOpt, build_options)
from sympy.polys.polyutils import (expr_from_dict, _dict_reorder,
                                  _parallel_dict_from_expr)
---> 31 from sympy.printing.defaults import DefaultPrinting
from sympy.utilities import public, subsets
from sympy.utilities.iterables import is_sequence

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/printing/__init__.py:11
from .mathml import mathml, print_mathml
from .python import python, print_python
---> 11 from .pycode import pycode
from .codeprinter import print_ccode, print_fcode
from .codeprinter import ccode, fcode, cxxcode, rust_code # noqa:F811

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/printing/pycode.py:11
from sympy.core.mod import Mod
from .precedence import precedence
---> 11 from .codeprinter import CodePrinter
_kw = {
   'and', 'as', 'assert', 'break', 'class', 'continue', 'def', 'del', 'elif',
   'else', 'except', 'finally', 'for', 'from', 'global', 'if', 'import', 'in',
   'is', 'lambda', 'not', 'or', 'pass', 'raise', 'return', 'try', 'while',
   'with', 'yield', 'None', 'False', 'nonlocal', 'True'
}
_known_functions = {
   'Abs': 'abs',
   'Min': 'min',
   'Max': 'max',
}

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/printing/codeprinter.py:13
from sympy.core.sorting import default_sort_key
from sympy.core.symbol import Symbol
---> 13 from sympy.functions.elementary.complexes import re
from sympy.printing.str import StrPrinter
from sympy.printing.precedence import precedence, PRECEDENCE

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/functions/__init__.py:50
from sympy.functions.special.elliptic_integrals import (elliptic_k,
       elliptic_f, elliptic_e, elliptic_pi)
from sympy.functions.special.beta_functions import beta, betainc, betainc_regularized
---> 50 from sympy.functions.special.mathieu_functions import (mathieus, mathieuc,
       mathieusprime, mathieucprime)
ln = log
__all__ = [
   'factorial', 'factorial2', 'rf', 'ff', 'binomial', 'RisingFactorial',
   'FallingFactorial', 'subfactorial',
   (...)    114     'mathieus', 'mathieuc', 'mathieusprime', 'mathieucprime',
]

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/functions/special/mathieu_functions.py:9
from sympy.functions.elementary.miscellaneous import sqrt
from sympy.functions.elementary.trigonometric import sin, cos
----> 9 class MathieuBase(DefinedFunction):
   """
   Abstract base class for Mathieu functions.

   This class is meant to reduce code duplication.

   """
   unbranched = True

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/core/basic.py:218, in Basic.__init_subclass__(cls)
def __init_subclass__(cls):
   # Initialize the default_assumptions FactKB and also any assumptions
   # property methods. This method will only be called for subclasses of
   # Basic but not for Basic itself so we call
   # _prepare_class_assumptions(Basic) below the class definition.
   super().__init_subclass__()
--> 218     _prepare_class_assumptions(cls)

File ~/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/sympy/core/assumptions.py:632, in _prepare_class_assumptions(cls)
defs = {}
for base in reversed(cls.__bases__):
--> 632     assumptions = getattr(base, '_explicit_class_assumptions', None)
   if assumptions is not None:
       defs.update(assumptions)

KeyboardInterrupt: 

s = HarmonicOscillator(n0=1.0)
s.plot()
assert abs(s.compute_dy_dt(s.copy())[...]).max() < 1e-14
assert np.allclose(s.get_energy(), -s.g * s.n0**2 * s.sigma * np.sqrt(np.pi / 2) / 2)

We can also do this in higher dimensions:

\[\begin{gather*} \psi_0(\vect{x}) = \sqrt{n_0}\exp\left(-\sum_{i}\frac{x_i^2}{2\sigma_i^2}\right), \qquad \frac{-\hbar^2}{2m}\nabla^2\psi_0 = \frac{\hbar^2}{2m}\sum_{i} \left(\frac{1}{\sigma_i^2} - \frac{x_i^2}{\sigma_i^4}\right) = \sum_{i} \left(\frac{\hbar^2}{2m\sigma_i^2} - \frac{m\omega_i^2x_i^2}{2}\right) \\ \qquad V(x) = \sum_i\left(\frac{m\omega_i^2x_i^2}{2} - \frac{\hbar^2}{2m\sigma_i^2}\right) - g\abs{\psi_0}^2, \qquad \sigma_{i}^2 = \frac{\hbar}{m\omega_i}. \end{gather*}\]

This potential has the non-linear piece subtracted so that the constructed solution is an exact solution with zero chemical potential (i.e. a stationary state) and energy:

\[\begin{gather*} E = \int_{-\infty}^{\infty}\d{x}\;\frac{-gn_0^2}{2}e^{-2x^2/\sigma^2} = -\sqrt{\frac{\pi}{2}}\frac{gn_0^2\sigma}{2} \end{gather*}\]

Travelling Waves#

The standard GPE admits a family of periodic traveling waves with an analytic form:

\[\begin{gather*} \psi(x, t) = e^{-\I \mu t /\hbar}\psi_v(x_v), \qquad x_v = x - vt \end{gather*}\]

where the wave is moving with speed \(v\). We base our solution on the presentation of Hoefer and El El:2016.

The GPE can be expressed in terms of \(\psi(x, t)\) or in terms of \(\psi_v(x_v) = e^{-\I\mu t/\hbar} \psi_v(x- vt)\):

\[\begin{gather*} \I\hbar\dot{\psi} = -\frac{\hbar^2}{2m}\pdiff[2]{\psi}{x} + g\abs{\psi}^2\psi,\qquad 0 = -\frac{\hbar^2\psi_v''}{2m} + v \I \hbar \psi_v' + (g\abs{\psi_v}^2 - \mu)\psi_v = \left(\frac{\op{p}^2}{2m} - v\op{p} + g\abs{\psi_v}^2 - \mu\right)\psi_v. \end{gather*}\]

Note: this is not the typical Galilean transformation for quantum mechanical systems which includes an extra factor of \(e^{\I (mvx +mv^2t/2)\hbar}\). This additional phase allows one to shift the momentum, completing the square of the kinetic energy, and removing the \(v\op{p}\) linear term at the expense of shifting the chemical potential. The present form, however, is more general, and works with arbitrary dispersion, so we maintain it.

Hoefer and El El:2016 set \(\hbar=m=g=1\) which defines the following dimensions:

\[\begin{gather*} m_\mathrm{unit} = m, \qquad l_\mathrm{unit} = \frac{mg}{\hbar^2}, \qquad t_\mathrm{unit} = \frac{m^3g^2}{\hbar^5}. \end{gather*}\]

Using the Madelung transform:

(6)#\[\begin{gather} \psi(x, t) = \sqrt{\rho(x, t)}e^{\I\phi(x, t)}, \qquad u(x, t) = \phi'(x, t),\\ \dot{\rho} + (\rho u)' = 0, \qquad (\rho u)_{,t} + \left(\rho u^2 + \frac{\rho^2}{2}\right)' = \frac{1}{4}\bigl(\rho\, (\log\rho)''\bigr)'. \tag{2.111} \end{gather}\]

The solution can be expressed as:

\[\begin{gather*} \DeclareMathOperator{\sn}{sn} \psi(x, t) = e^{-\I\mu t/\hbar}\psi_v(x_v), \qquad \psi_v(x_v) = \sqrt{\rho_v(x_v)}e^{\I\phi_v(x_v)}, \qquad u(x, t) = \phi_v'(x, t),\\ m = al^2, \qquad \rho_v(x_v) = \rho_{\min} + a \sn^2\bigl(\frac{x_v}{l};m\bigr), \qquad x_v = x - vt, \qquad u(x, t) = v - \frac{C}{\rho(x, t)}, \\ a = \rho_{\max} - \rho_{\min}, \qquad \mu = \rho_{\min} - \frac{v^2}{2} + \frac{\rho_\max + a/m}{2}\qquad C = \frac{\sqrt{\rho_{\min}\rho_{\max}(1+\rho_{\min}l^2)}}{l},\qquad L = 2lK(m), \qquad Q = \frac{2\pi}{L}. \end{gather*}\]

(Note: \(k=\sqrt{m}\) for \(\sn(z,k)\) in some CASs. We use \(\sn(z;m)\) and \(K(m)\) here)

The full solution (with proper coefficients) is thus:

\[\begin{gather*} \psi(x, t) = e^{-\I\mu t/\hbar}\psi_v(x_v), \qquad \psi_v(x_v) = \sqrt{\rho_v(x_v)}e^{\I\phi_v(x_v)}, \qquad u(x_v) = \frac{\hbar}{m}\phi_v'(x_v) = v - \frac{C}{\rho_v(x_v)},\\ \rho_v(x_v) = \rho_{\min} + a\sn^2\bigl(\frac{x_v}{l};\tfrac{mg}{\hbar^2}al^2\bigr), \qquad x_v = x - vt, \qquad a = \rho_{\max} - \rho_{\min}, \\ \mu = g\rho_\min - \frac{mv^2}{2} + g\frac{\rho_\max +a/m}{2}\qquad C = \frac{\sqrt{\rho_{\min}\rho_{\max}(\hbar^2+mgl^2\rho_{\min})}}{ml}. \end{gather*}\]

\[\begin{gather*} m = al^2, \qquad l^{-1} = \frac{2K(m)}{L} \end{gather*}\]

%pylab inline --no-import-all
from gpe.imports import *
import gpe.exact_solutions

reload(gpe.exact_solutions)
from gpe.exact_solutions import TravellingWaves

# Stationary wave
args = dict(n0=0.9, n1=1.0, Lx=10.0, Nx=16)
s0 = TravellingWaves(v_p=0, **args)

# Travelling wave in the lab frame
s1 = TravellingWaves(v_p=0.6047197603, v_x=0, **args)

# Travelling wave in a co-moving frame.
s2 = TravellingWaves(v_p=0.6047197603, **args)

for ss in evolves([s0, s1, s2], t_max=50.0):
    plt.clf()
    for s in ss:
        s.plot()

# s.twist_x
# s.mu, s.get_mu().real, s._twist

As a simple test, we consider the phonon limit \(a = \epsilon \ll 1\). In this limit, we have

\[\begin{gather*} K(m) \approx \frac{\pi}{2}\left(1 + \frac{m}{4}\right) + \order(m^2), \qquad m = \frac{L^2}{\pi^2}\epsilon, \qquad \sn(\theta, m) = \sin(\theta) + \order(m), \qquad l = \frac{L}{\pi}\left(1 - \frac{L^2}{4\pi^2}\epsilon\right) + \order(\epsilon^2)\\ \end{gather*}\]

\[\begin{gather*} \frac{L = \approx l\pi\left(1 + \frac{m}{4}\right) \end{gather*}\]

consider a \(v=0\) solution with \(\rho_\min = \rho_\max = 1\):

\[\begin{gather*} \rho = 1, \qquad a = m = 0, \qquad \frac{L}{l} = \pi, \qquad u(x_v) = \frac{C}{\rho}, \qquad \theta(x_v) = \frac{C}{\rho}x_v \end{gather*}\]

from gpe.exact_solutions import K

s._C

s = TravellingWaves(n0=1.0, n1=1.0, Lx=np.pi, Nx=64, v_p=0.0)
s.plot()
s.twist_x
s._mu, s.get_mu().real, s._twist

e = EvolverABM(s, dt=0.2 * s.t_scale)
with NoInterrupt() as interrupted:
    while not interrupted:
        e.evolve(1000)
        plt.clf()
        e.y.plot()
        display(plt.gcf())
        plt.close("all")
        clear_output(wait=True)

x = s.xyz[0]
c = np.sqrt(s.g * s.n1 / s.m)
v = 0.5 * c
u = np.sqrt(c**2 - v**2)
n0 = s.n1 * (v**2 / c**2)
l = s.hbar / s.m / u
psi_soliton = v / c * 1j + u / c * np.tanh(x / l)
twist = np.angle(psi_soliton)[-1] - np.angle(psi_soliton)[0]
s1 = TravellingWaves(n0=n0, n1=s.n1, Lx=s.basis.Lx, Nx=s.basis.Nx, v_p=0, twist=twist)
s1[...] = psi_soliton
s1.plot()

e = EvolverABM(s, dt=0.2 * s.t_scale)
e1 = EvolverABM(s1, dt=0.2 * s1.t_scale)
pe = None
with NoInterrupt() as interrupted:
    while not interrupted:
        e.evolve(1000)
        e1.evolve(1000)
        plt.clf()
        e.y.plot()
        e1.y.plot()
        display(plt.gcf())
        plt.close("all")
        clear_output(wait=True)

Bright Soliton#

Here we demonstrate the analytic bright-soliton for a GPE with attractive interactions. This object moves at a specified speed \(v\), and we consider the solution in three frames: comoving (the soliton is stationary), lab frame (the soliton makes a single oscillation in time \(T=L/v\)) and a frame moving backwards at the same speed (the soliton crosses the box twice in time \(T\)).

This forms the basis of gpe/tests/test_bec.py::test_StateTwist_x_v_x and verifies that the frame velocity argument of StateTwist_x works. Because the density vanishes at the boundaries, this does not test the twist.

\[\begin{gather*} \I\hbar\dot{\psi}(x, t) = \frac{-\hbar^2\psi''(x, t)}{2m} + g\abs{\psi}^2\psi,\\ \psi(x, t) = \frac{\sqrt{n_0}}{\cosh\bigl(\eta (x-vt)\bigr)} \exp\left\{\frac{\left(\mu + \tfrac{mv^2}{2}\right)t - vx}{\I\hbar}\right\}, \qquad \mu = \frac{-\hbar^2\eta^2}{2m}, \\ gn_0 = 2\mu = \frac{-\hbar^2\eta^2}{m}, \qquad \sigma = \frac{1}{\eta}. \end{gather*}\]

import sympy
from sympy import Eq, var, I, exp, cosh, sqrt, Abs

x, t, g = var(["x", "t", "g"], real=True)
eta, n0, mu, m, hbar = var(["eta", "n_0", "mu", "m", "hbar"], positive=True)
mu = -(hbar**2) * eta**2 / 2 / m
g = 2 * mu / n0
psi = exp(mu * t / I / hbar) * sqrt(n_0) / cosh(eta * x)
n = Abs(psi) ** 2
display(Eq(sympy.S("psi"), psi))
Hpsi = -(hbar**2) * psi.diff(x, x) / 2 / m + g * n * psi
display(Eq(sympy.S("H*psi"), Hpsi))
(Hpsi - I * hbar * psi.diff(t)).simplify()

%pylab inline --no-import-all
from gpe.imports import *
import gpe.exact_solutions

reload(gpe.exact_solutions)
from gpe.exact_solutions import BrightSoliton

Nx = 32
Lx = 10.0
v = 2.0

args = dict(Nx=Nx, Lx=Lx, v=v, sigma=1)
s0 = BrightSoliton(v_x=0, **args)
s1 = BrightSoliton(v_x=-v, **args)
s2 = BrightSoliton(v_x=v, **args)

for ss in evolves([s0, s1, s2], t_max=Lx / v):
    plt.clf()
    for s in ss:
        s.plot()

1D GPE (NLSEQ)#

Here we consider the exact solutions to the 1D GPE via the inverse scattering method (ISM). These follow from a so-called Lax representation, which can be expressed as follows. Consider two linear operators \(\op{L}(\lambda)\) and

\[\begin{gather*} \op{M}(\lambda) = \I\diff{}{t} + V(x, t, \lambda) \end{gather*}\]

that commute:

\[\begin{gather*} [\op{L}(\lambda), \op{M}(\lambda)] = 0. \end{gather*}\]

\[\begin{gather*} \mat{L} = \I\mat{1}\pdiff{}{x} + \begin{pmatrix} 0 & u^*\\ u & 0 \end{pmatrix},\\ \mat{M} = -\mat{1}\pdiff[2]{}{x} + \begin{pmatrix} \abs{u}^2 & \I u_{,x}^*\\ -\I u_{,x} & -\abs{u}^2 \end{pmatrix}. \end{gather*}\]

differential equation

\[\begin{gather*} \I \dot{\mat{L}} = [\mat{L}, \mat{A}] \end{gather*}\]

We start by expressing the problem as

\[\begin{gather*} \I \dot{u} + u'' + 2 \abs{u}^2 u = 0, \qquad u(x, t=0) = u_0(x). \end{gather*}\]

Define

\[\begin{gather*} \mat{q}_0(x) = \begin{pmatrix} 0 & u_0(x)\\ u_0^*(x) & 0 \end{pmatrix}. \end{gather*}\]