# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Julia 0.5.2
#     language: julia
#     name: julia-0.5
# ---

# ### Background  
# 1. C-field theory
#  * Projected fields
#  * Projected Gross-Pitaevskii Equation
# 2. Spectral-Galerkin method
#  * $2M$ Rule for quadrature integration
#  * Integrating general nonlinearities

# ### C-field definition
# \begin{equation}
# \psi(x,t)=\sum_n c_n(t)\phi_n(x),
# \end{equation}
# where the modes $\phi_n(x)$ are a complete and orthonormal set:
# \begin{equation}
# \int dx\;\phi^*_n(x)\phi_m(x)=\delta_{nm},\\
# \mathbb{1}=\sum_n|\phi_n\rangle\langle\phi_n|.
# \end{equation}
# Completeness then gives a construction for the Dirac delta:
# \begin{equation}
# \delta(x-x')=\sum_n\phi_n(x)\phi_n^*(x').
# \end{equation}

# ### Defining the C-region
# To define the C-region, we introduce the projection operator
# \begin{equation}
# \hat{\cal P}=\bar{\sum_n}|\phi_n\rangle\langle \phi_n|
# \end{equation}
#
# where the bar notation refers to the cutoff defining the $C$-region. The quantum number is restricted to $n\leq M-1$, giving a total of $M$ modes in the 1D $C$-region, indexed by quantum numbers $n\in C=\{ 0,1,\dots, M-1\}$.
#
# The Dirac delta function for the $C$-region is 
# \begin{equation}
# \langle x|\hat{\cal P}|x'\rangle=\bar{\sum_n}\phi_n(x)\phi_n^*(x')\equiv\delta(x,x').
# \end{equation}

# ### Projected $\delta$-function
# The projected delta function $\delta(x,x')$ plays a fundamental role in  c-field theory. It forms the kernel of the projector acting on functions in  position space:
# \begin{equation}
# {\cal P}f(x)\equiv \int dx'\;\delta(x,x')f(x')
# \end{equation}
# For any function residing in $C$, projection is the identity. Equivalently, the action of $\delta(x,x')$ is equivalent to the defining property of $\delta(x-x')$:
# \begin{equation}
# \int dx'\;\delta(x,x')f(x')=f(x).
# \end{equation}
# However, when considering $\delta(x,x')$, the order of arguments is significant.

# ## Projected  Gross-Pitaevskii equation
# ### Introduction to PGPE
# Given a basis as above, one start from the GPE
# \begin{equation}
# i\hbar\frac{\partial \psi}{\partial t}=\left(-\frac{\hbar^2\nabla^2}{2m}+V(x)+g|\psi|^2\right)\psi
# \end{equation}
# and substitute 
# \begin{equation}
# \psi=\bar{\sum_n}c_n(t)\phi_n(x),
# \end{equation}
# and project onto $\phi_n(x)$.

# ## Spectral basis
# We express the equation in the spectral basis 
# \begin{equation}
# (\hat T+\hat V)\phi_n(x)=\epsilon_n\phi_n(x).
# \end{equation}
# We arrive at the equation
# \begin{equation}
# i\hbar\frac{\partial c_n}{\partial t}=\epsilon_nc_n+g\int dx\;\phi_n^*(x)|\psi(x)|^2\psi(x)
# \end{equation}
# In the harmonic oscillator basis, our modes take the form 
# \begin{equation}
# \phi_n(x)=e^{-x^2/2}P_n(x)
# \end{equation}
# where $P_n(x)$ contains the polynomial part of degree $n$ and normalization factors. 

# There are two obvious ways to express the integral.
# #### 1. Direct substitution of basis decomposition
# \begin{equation}
# I_n\equiv \int dx\;\phi_n^*(x)|\psi(x)|^2\psi(x)
# \end{equation}
# Firstly, we can substitute the wavefunction:
# \begin{equation}
# I_n=\bar{\sum}_{mpq}c_m^*c_pc_q\int dx\;\phi_n^*(x)\phi_m^*(x)\phi_p(x)\phi_q(x),
# \end{equation}
# a triple summation over a set of known matrix elements for each mode $n$. For each $n$ this triple summation  requires $O(M^3)$ arithmetic operations. The total cost of evaluating the PGPE right hand side via summation is $O(M^4)$. This approach is numerically expensive, and scales poorly with $M$.

# ####  2. Quadrature method
# This computational effort can be reduced via an alternate approach using Gaussian quadrature. 
#
# >#### Philosophy:   
# We will transform onto a *particular* (in general non-uniform) spatial grid that will guarantee that the non-linear term is projected faithfully onto $\phi_n(x)$, *for all modes in the C-field*.
#
# For the integral at hand, the highest order term can be found by considering the cutoff mode for each wavefunction $\psi(x)$, and can be written as
# \begin{equation}
# I_{M-1}=\int dx\;e^{-2x^2}R_{4(M-1)}(x),
# \end{equation}
# where $R$ expresses the polynomial part.

# ### The $2M$ rule
# As we will see below, for this $C$ region consisting of $M$ modes, the nonlinear integral may be evaluated exactly by mapping the coefficients $c_n$ to a spatial grid of $2M$ points, corresponding to a Gaussian quadrature rule of order  $2M$. Thus we can state our general rule, for a cubic nonlinearity of Gross-Pitaevskii form:
#
# >#### The 2M rule:
# A field represented by $M$ modes may be evolved at working precision using a particular grid of $2M$ spatial points.
#
# The Fourier spectral method is a special case for which the wieght function becomes unity, and for which the $2M$ rule achieves spatial sampling at the *Nyquist frequency*.

# ### Gauss-Hermite quadrature
# For any integral of the form
# \begin{equation}
# I=\int_{-\infty}^{\infty} dx \;e^{-x^2}Q_{2N-1}(x),
# \end{equation}
# with $Q_{2N-1}(x)$ restricted to maximum polynomial degree $2N-1$
# \begin{equation}
# Q_{2N-1}(x)=a_0+a_1x+\dots+a_{2N-1}x^{2N-1},
# \end{equation}
# there exists a quadrature rule  of order $N$, involving $N$ roots $x_k$, and $N$ weights $w_k$, $k\in 1,\dots, N$, that will evaluated all such integrals exactly. 
#
# From a linear algebra point of view this result is not so surprising: The Gaussian-weighted integral of an arbitrary polynomial of degree $2N-1$ involves $2N$ uknown coefficients, and $2N$ known integrals. 

# Given $2N$ free parameters, here given by the roots and weights of the quadrature rule, the $2N$ unknowns may be solved for exactly, by solving a system of simultaneous linear equations. 
#
# The integral may be evaluated "exactly" as
# \begin{equation}
# I=\sum_k w_kQ_{2N-1}(x_k),
# \end{equation}
# where the meaning of  "exact" is that it is accurate to machine precision; typically the accuracy will be of order 1e-16. 

# ### Projected time evolution: work step
# #### Algorithm
# >1. cast the integral in the above form $I$, 
# 2. Transform to the quadrature grid $c_n\to \psi(x_k)$
# 2. Evaluate $Q_{2N-1}(x_k)$ for the resulting polynomial part, and 
# 3. Evaluate the sum weighted by $w_k$, giving the projection of $|\psi|^2\psi\to c_n'$
#
# Step 1. is done by the change of variables $x=y/\sqrt{2}$

# we find
# \begin{align}
# I_{M-1}&=\int_{-\infty}^{\infty} dy\;e^{-y^2}\left[\frac{R_{4(M-1)}(y/\sqrt{2})}{\sqrt{2}}\right],
# \end{align}
# where the terms in parentheses correspond to $Q_{2N-1}(x)$ defined above. 
#
# According to the rules of Gaussian quadrature, a rule of order $N$ will integrate all polynomials up to and including degree $2N-1$. A rule of order $2M$ will thus evaluate all terms up to order $4M-1$. This is the lowest order rule that will guarantee exactness.
# Hence
# >**All projected integrals in the PGPE for $M$ modes are exact within a Gaussian quadrature of order $2M$**.

# A rule of order $n$ will integrate polynomials up to and including degree $2n-1$. Thus a rule of order $JM$ will integrate terms up to order $2JM-1$, giving numerically exact evalution of the nonlinear term for all modes in $C$. The variable transformation is 
# $$ x=\frac{y}{\sqrt{J}}$$
# giving the integral 
# $$I_{M-1}=\int_{-\infty}^{\infty} dy\;e^{-y^2} \frac{1}{\sqrt{J}}R_{2J(M-1)}(y/\sqrt{J}),$$

# All powers of $\psi(x)$, i.e. both odd and even, are accommodated by putting $2J\to K$, to give
# $$
# I_{M-1}\equiv \int dx\;\phi_n^*(x)\psi(x)^{K-2}\psi(x),
# $$
# to give the form
# $$
# I_{M-1}=\int_{-\infty}^{\infty} dx \;e^{-Kx^2/2}R_{K(M-1)}(x).
# $$
# A rule of order $KM/2$ will integrate all terms up to order $KM-1>K(M-1)$. The variable transformation is 
# $$x=\sqrt{\frac{2}{K}}y,$$ 
# giving the integral
# $$I_{M-1}=\int_{-\infty}^{\infty} dy\;e^{-y^2} \sqrt{\frac{2}{K}}R_{K(M-1)}\left(\sqrt{\frac{2}{K}}y\right) ,$$

# ### General Anisotropic form
# The most general formulation we require must accommodate different oscillator fequencies in different spatial directions:
# $$
# I_{M-1}\equiv \int dx\;\phi_n^*(x)\psi(x)^{K-2}\psi(x),
# $$
# to give the form
# $$
# I_{M-1}=\int_{-\infty}^{\infty} dx \;e^{-K\omega x^2/2}R_{K(M-1)}(x).
# $$
# A rule of order $KM/2$ will integrate all terms with polynomial degree up to $KM-1>K(M-1)$. The variable transformation is 
# $$x=\sqrt{\frac{2}{K\omega}}y,$$ 
# giving the integral
# $$I_{M-1}=\int_{-\infty}^{\infty} dy\;e^{-y^2} \sqrt{\frac{2}{K\omega}}R_{K(M-1)}\left(\sqrt{\frac{2}{K\omega}}y\right).$$
#
# In the package `ProjectedGPE.jl`, the function that sets up the roots `x`, weights `w`, and transformation matrix `T` to evaluate integrals of this form is 
# `nfieldtrans.jl`

# You should be able to run this example to compute wavefunction norm:

using ProjectedGPE;

M = 20
x,w,T = nfieldtrans("Hermite",M,2) #2 field product in the integral to be evaluated
c = randn(M) + im*randn(M);
ψ = T*c
N = sum(w.*abs(ψ).^2)

# and compare with the direct summation of coefficients:

Nsum = sum(abs(c).^2)

abs(N-Nsum)

# Integrals of higher order field products are computed with comparable accuracy.



StarcoderdataJupyterScriptsDedupFiltered