notebook.community

Edit and run

Table of Contents

  1. Rabi Oscillations
    1. Simulating the Full Hamiltonian
    2. With Rotating Wave Approximation
  2. Coherent State in a Harmonic Oscillator
  3. Jaynes-Cummings Revival
    1. Definite Photon State
    2. Coherent State

The only addition to this notebook compared to the original is the DEBUG=True setting that prints all of the generated Cython code.





In [1]:



    


from cutiepy import *
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
cutiepy.codegen.DEBUG = True

Rabi Oscillations

Simulating the Full Hamiltonian

$\hat{H} = \hat{H}_0 + \Omega \sin((\omega_0+\Delta)t) \hat{\sigma}_x$

$\hat{H}_0 = \frac{\omega_0}{2}\hat{\sigma}_z$





In [2]:



    


initial_state = basis(2, 0)
initial_state



    


















    
Out[2]:







$$\text{Ket }{| {{0}}_{{\tiny N\normalsize 2}} \rangle} \text{ on the space }\mathbb{C}^{2}\text{ with numerical content: }$$
$$\begin{equation*}\left(\begin{array}{*{11}c}1.0\\0.0\\\end{array}\right)\end{equation*}$$
















In [3]:



    


ω0 = 1
Δ = 0.002
Ω = 0.005
ts = 2*np.pi/Ω*np.linspace(0,1,40)
H = ω0/2 * sigmaz() + Ω * sigmax() * sin((ω0+Δ)*t)
H



    


















    
Out[3]:







$${\left( {{0.005}\tiny\times\normalsize{\operatorname{sin}\left( {{{1.002}\tiny\times\normalsize{t}}} \right)}\tiny\times\normalsize{\hat{σ}_{x}}}+{{0.5}\tiny\times\normalsize{\hat{σ}_{z}}} \right)}$$
















In [4]:



    


res = sesolve(H, initial_state, ts)



    


















    






Generating cython code...
Compiling cython code...










    











  
  
  






#cython: boundscheck=False, wraparound=False, initializedcheck=False, cdivision=True
cimport numpy as np
import numpy as np
from libc.math cimport sin, cos, exp, sinh, cosh, tan, tanh, log
from scipy.linalg.cython_blas cimport zgemv, zgemm

cdef int iONE=1, iZERO=0
cdef double complex zONE=1, zZERO=0, zNHALF=-0.5

# Declaration of global variables for
# - intermediate results
# - predefined constants
cdef double D = 0
cdef double B = 0
cdef double complex [:, :] F = np.empty((2, 2), dtype=np.complex)
cdef double complex [:, :] G = np.empty((2, 2), dtype=np.complex)
cdef double complex [:, :] I = np.empty((2, 1), dtype=np.complex)
cdef double complex [:, :] E # initialized in setup
cdef double complex [:, :] A # initialized in setup

# The generated cython code
cdef generated_function(
    # input arguments
    double C,
    double complex [:, :] H,
    # output argument
    double complex [:, :] J
    ):
    # indices and intermediate values for various matrix operations
    cdef int i, j, k, ii, jj, kk
    cdef double complex c
    # intermediate results
    global D, B, F, G, I
    # predefined constants
    global E, A
    # evaluating the function
    D = 1.002*C
    B = sin(D)
    # F = B*E
    ii = E.shape[0]
    jj = E.shape[1]
    for i in range(ii):
        for j in range(jj):
            F[i,j] = B*E[i,j]
    # G = A+F
    ii = A.shape[0]
    jj = A.shape[1]
    for i in range(ii):
        for j in range(jj):
            G[i,j] = A[i,j]+F[i,j]
    # I = G.H
    ii = G.shape[0]
    jj = G.shape[1]
    kk = H.shape[1]
    # Commented out col-major order (BLAS default):
   #zgemm('N', 'N', &ii, &kk, &jj, &zONE, &G[0,0], &ii, &H[0,0], &jj, &zZERO, &I[0,0], &ii)
    zgemm('N', 'N', &kk, &ii, &jj, &zONE, &H[0,0], &kk, &G[0,0], &jj, &zZERO, &I[0,0], &kk)
    # J = -1j*I
    ii = I.shape[0]
    jj = I.shape[1]
    for i in range(ii):
        for j in range(jj):
            J[i,j] = -1j*I[i,j]

# Function to deliver the constants from python to cython
cpdef setup_generated_function(
    double complex [:, :] _E,
    double complex [:, :] _A
    ):
    global E, A
    E = _E
    A = _A

ctypedef void rhs_t(double t, double *y, double *ydot)
cdef extern from "cvode_simple_interface.h":
    struct cvsi_instance:
        pass
    cvsi_instance *cvsi_setup(rhs_t *rhs,
                              double *y0, int neq,
                              double t0,
                              long mxsteps, double reltol, double abstol)
    int cvsi_step(cvsi_instance *instance, double t)
    void cvsi_destroy(cvsi_instance *instance)

cdef inline void RHS(double t, double *y, double *ydot):
    generated_function(t,
                       <double complex [:2, :1]> <double complex *> y,
                       <double complex [:2, :1]> <double complex *> ydot)
    # TODO The casts above need to allocate new memview structs each time
    # It is noticeable on 2x2 matrices (10%), not noticeable otherwise

cpdef list pythonsolve(
        np.ndarray[np.double_t, ndim=1] ts,
        np.ndarray[np.complex_t, ndim=2] y0,
        long mxsteps, double rtol, double atol,
        progressbar):
    if ts[0] == 0:
        ts = ts[1:]
        res = [y0]
    else:
        res = []
    cdef np.ndarray[np.complex_t, ndim=2, mode="c"] _y0 = np.copy(y0)
    cdef int neq = np.size(_y0)*2
    cdef cvsi_instance *instance = cvsi_setup(RHS, <double *> &_y0[0,0], neq,
                                              0, mxsteps, rtol, atol)
    for t in ts:
        cvsi_step(instance, t)
        res.append(np.copy(_y0))
        progressbar.step()
    cvsi_destroy(instance)

    progressbar.stop()
    return res














    






Running cython code...
Starting at 10/20 13:52:34.
Finishing at 10/20 13:52:34.
Total time: 0 seconds.
Formatting the output...

















In [5]:



    


σz_expect = expect(sigmaz(), res)



    











In [6]:



    


plt.plot(ts*Ω/np.pi, σz_expect, 'r.', label='numerical result')
Ωp = (Ω**2+Δ**2)**0.5
plt.plot(ts*Ω/np.pi, 1-(Ω/Ωp)**2*2*np.sin(Ωp*ts/2)**2, 'b-',
         label=r'$1-2(\Omega^\prime/\Omega)^2\sin^2(\Omega^\prime t/2)$')
plt.title(r'$\langle\sigma_z\rangle$-vs-$t\Omega/\pi$ at '
          r'$\Delta/\Omega=%.2f$, $\omega_0/\Omega=%.2f$'%(Δ/Ω, ω0/Ω))
plt.ylim(-1,1)
plt.legend(loc=3);

With Rotating Wave Approximation

$\hat{H}^\prime = e^{i\hat{H}_0 t}\hat{H} e^{-i\hat{H}_0 t} \approx \frac{\Delta}{2} \hat{\sigma}_z + \frac{\Omega}{2} \hat{\sigma}_x$





In [7]:



    


Hp = Δ/2 * sigmaz() + Ω/2 * sigmax()
Hp



    


















    
Out[7]:







$${\left( {{0.001}\tiny\times\normalsize{\hat{σ}_{z}}}+{{0.0025}\tiny\times\normalsize{\hat{σ}_{x}}} \right)}$$
















In [8]:



    


res = sesolve(Hp, initial_state, ts)



    


















    






Generating cython code...
Compiling cython code...










    











  
  
  






#cython: boundscheck=False, wraparound=False, initializedcheck=False, cdivision=True
cimport numpy as np
import numpy as np
from libc.math cimport sin, cos, exp, sinh, cosh, tan, tanh, log
from scipy.linalg.cython_blas cimport zgemv, zgemm

cdef int iONE=1, iZERO=0
cdef double complex zONE=1, zZERO=0, zNHALF=-0.5

# Declaration of global variables for
# - intermediate results
# - predefined constants
cdef double complex [:, :] C = np.empty((2, 1), dtype=np.complex)
cdef double complex [:, :] A # initialized in setup

# The generated cython code
cdef generated_function(
    # input arguments
    double E,
    double complex [:, :] B,
    # output argument
    double complex [:, :] D
    ):
    # indices and intermediate values for various matrix operations
    cdef int i, j, k, ii, jj, kk
    cdef double complex c
    # intermediate results
    global C
    # predefined constants
    global A
    # evaluating the function
    # C = A.B
    ii = A.shape[0]
    jj = A.shape[1]
    kk = B.shape[1]
    # Commented out col-major order (BLAS default):
   #zgemm('N', 'N', &ii, &kk, &jj, &zONE, &A[0,0], &ii, &B[0,0], &jj, &zZERO, &C[0,0], &ii)
    zgemm('N', 'N', &kk, &ii, &jj, &zONE, &B[0,0], &kk, &A[0,0], &jj, &zZERO, &C[0,0], &kk)
    # D = -1j*C
    ii = C.shape[0]
    jj = C.shape[1]
    for i in range(ii):
        for j in range(jj):
            D[i,j] = -1j*C[i,j]

# Function to deliver the constants from python to cython
cpdef setup_generated_function(
    double complex [:, :] _A
    ):
    global A
    A = _A

ctypedef void rhs_t(double t, double *y, double *ydot)
cdef extern from "cvode_simple_interface.h":
    struct cvsi_instance:
        pass
    cvsi_instance *cvsi_setup(rhs_t *rhs,
                              double *y0, int neq,
                              double t0,
                              long mxsteps, double reltol, double abstol)
    int cvsi_step(cvsi_instance *instance, double t)
    void cvsi_destroy(cvsi_instance *instance)

cdef inline void RHS(double t, double *y, double *ydot):
    generated_function(t,
                       <double complex [:2, :1]> <double complex *> y,
                       <double complex [:2, :1]> <double complex *> ydot)
    # TODO The casts above need to allocate new memview structs each time
    # It is noticeable on 2x2 matrices (10%), not noticeable otherwise

cpdef list pythonsolve(
        np.ndarray[np.double_t, ndim=1] ts,
        np.ndarray[np.complex_t, ndim=2] y0,
        long mxsteps, double rtol, double atol,
        progressbar):
    if ts[0] == 0:
        ts = ts[1:]
        res = [y0]
    else:
        res = []
    cdef np.ndarray[np.complex_t, ndim=2, mode="c"] _y0 = np.copy(y0)
    cdef int neq = np.size(_y0)*2
    cdef cvsi_instance *instance = cvsi_setup(RHS, <double *> &_y0[0,0], neq,
                                              0, mxsteps, rtol, atol)
    for t in ts:
        cvsi_step(instance, t)
        res.append(np.copy(_y0))
        progressbar.step()
    cvsi_destroy(instance)

    progressbar.stop()
    return res














    






Running cython code...
Starting at 10/20 13:52:58.
Finishing at 10/20 13:52:58.
Total time: 0 seconds.
Formatting the output...

















In [9]:



    


σz_expect = expect(sigmaz(), res)



    











In [10]:



    


plt.plot(ts*Ω/np.pi, σz_expect, 'r.', label='numerical result')
Ωp = (Ω**2+Δ**2)**0.5
plt.plot(ts*Ω/np.pi, 1-(Ω/Ωp)**2*2*np.sin(Ωp*ts/2)**2, 'b-',
         label=r'$1-2(\Omega^\prime/\Omega)^2\sin^2(\Omega^\prime t/2)$')
plt.title(r'$\langle\sigma_z\rangle$-vs-$t\Omega/\pi$ at '
          r'$\Delta/\Omega=%.2f$ in RWA'%(Δ/Ω))
plt.ylim(-1,1)
plt.legend(loc=3);

Coherent State in a Harmonic Oscillator

$|\alpha\rangle$ evolving under $\hat{H} = \hat{n}$





In [11]:



    


N_cutoff = 40
α = 2.5
initial_state = coherent(N_cutoff, α)
initial_state



    


















    
Out[11]:







$$\text{Ket }{| {{\tiny\alpha\normalsize 2.50}}_{{\tiny N\normalsize 40}} \rangle} \text{ on the space }\mathbb{C}^{40}\text{ with numerical content: }$$
$$\begin{equation*}\left(\begin{array}{*{11}c}0.044\\0.110\\0.194\\0.280\\0.350\\\vdots\\3.661\times10^{-08}\\1.525\times10^{-08}\\6.270\times10^{-09}\\2.543\times10^{-09}\\1.018\times10^{-09}\\\end{array}\right)\end{equation*}$$
















In [12]:



    


H = num(N_cutoff)
H



    


















    
Out[12]:







$$\text{Operator }{\hat{{n}}_{{40}}} \text{ on the space }\mathbb{C}^{40}\text{ with numerical content: }$$
$$\begin{equation*}\left(\begin{array}{*{11}c}0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 1.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 2.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 3.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 4.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\\vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \vdots\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 35.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 36.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 37.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 38.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 39.0\\\end{array}\right)\end{equation*}$$
















In [13]:



    


ts = 2*np.pi*np.linspace(0,1,41)
res = sesolve(H, initial_state, ts)
a = destroy(N_cutoff)
a_expect = expect(a, res, keep_complex=True)



    


















    






Generating cython code...
Compiling cython code...










    











  
  
  






#cython: boundscheck=False, wraparound=False, initializedcheck=False, cdivision=True
cimport numpy as np
import numpy as np
from libc.math cimport sin, cos, exp, sinh, cosh, tan, tanh, log
from scipy.linalg.cython_blas cimport zgemv, zgemm

cdef int iONE=1, iZERO=0
cdef double complex zONE=1, zZERO=0, zNHALF=-0.5

# Declaration of global variables for
# - intermediate results
# - predefined constants
cdef double complex [:, :] C = np.empty((40, 1), dtype=np.complex)
cdef double complex [:] A # initialized in setup
cdef int [:] A_pointers # initialized in setup
cdef int [:] A_indices # initialized in setup

# The generated cython code
cdef generated_function(
    # input arguments
    double E,
    double complex [:, :] B,
    # output argument
    double complex [:, :] D
    ):
    # indices and intermediate values for various matrix operations
    cdef int i, j, k, ii, jj, kk
    cdef double complex c
    # intermediate results
    global C
    # predefined constants
    global A, A_pointers, A_indices
    # evaluating the function
    # C = A.B
    ii = A_pointers.shape[0]-1
    for i in range(ii):
        jj = A_pointers[i]
        kk = A_pointers[i+1]
        c = 0
        for j in range(jj,kk):
            k = A_indices[j]
            c += A[j]*B[k,0]
        C[i,0] = c
    
    # D = -1j*C
    ii = C.shape[0]
    jj = C.shape[1]
    for i in range(ii):
        for j in range(jj):
            D[i,j] = -1j*C[i,j]

# Function to deliver the constants from python to cython
cpdef setup_generated_function(
    double complex [:] _A,
    int [:] _A_pointers,
    int [:] _A_indices
    ):
    global A, A_pointers, A_indices
    A = _A
    A_pointers = _A_pointers
    A_indices = _A_indices

ctypedef void rhs_t(double t, double *y, double *ydot)
cdef extern from "cvode_simple_interface.h":
    struct cvsi_instance:
        pass
    cvsi_instance *cvsi_setup(rhs_t *rhs,
                              double *y0, int neq,
                              double t0,
                              long mxsteps, double reltol, double abstol)
    int cvsi_step(cvsi_instance *instance, double t)
    void cvsi_destroy(cvsi_instance *instance)

cdef inline void RHS(double t, double *y, double *ydot):
    generated_function(t,
                       <double complex [:40, :1]> <double complex *> y,
                       <double complex [:40, :1]> <double complex *> ydot)
    # TODO The casts above need to allocate new memview structs each time
    # It is noticeable on 2x2 matrices (10%), not noticeable otherwise

cpdef list pythonsolve(
        np.ndarray[np.double_t, ndim=1] ts,
        np.ndarray[np.complex_t, ndim=2] y0,
        long mxsteps, double rtol, double atol,
        progressbar):
    if ts[0] == 0:
        ts = ts[1:]
        res = [y0]
    else:
        res = []
    cdef np.ndarray[np.complex_t, ndim=2, mode="c"] _y0 = np.copy(y0)
    cdef int neq = np.size(_y0)*2
    cdef cvsi_instance *instance = cvsi_setup(RHS, <double *> &_y0[0,0], neq,
                                              0, mxsteps, rtol, atol)
    for t in ts:
        cvsi_step(instance, t)
        res.append(np.copy(_y0))
        progressbar.step()
    cvsi_destroy(instance)

    progressbar.stop()
    return res














    






Running cython code...
Starting at 10/20 13:53:23.
Finishing at 10/20 13:53:23.
Total time: 0 seconds.
Formatting the output...

















In [14]:



    


plt.figure(figsize=(4,4))
plt.plot(np.real(a_expect), np.imag(a_expect), 'b-')
for t, alpha in list(zip(ts,a_expect))[:40:4]:
    plt.plot(np.real(alpha), np.imag(alpha), 'r.')
    plt.text(np.real(alpha), np.imag(alpha), r'$t=%.1f\pi$'%(t/np.pi), fontsize=14)
plt.title(r'$\langle\hat{a}\rangle$-vs-$t$')
plt.ylabel(r'$\mathcal{I}(\alpha)$')
plt.xlabel(r'$\mathcal{R}(\alpha)$');

Jaynes-Cummings Revival

$\hat{H} = \hat{H}_0 + \hat{H}^\prime$

$\hat{H}_0 = \omega \hat{n} + \omega \frac{1}{2} \hat{\sigma}_z$

$\hat{H}^\prime = g (\hat{\sigma}_+\hat{a} + \hat{\sigma}_-\hat{a}^\dagger)$





In [15]:



    


ω = 1
g = 0.1
ts = np.pi/g*np.linspace(0,1,150)
N_cutoff = 40
H0 = ω*(tensor(num(N_cutoff), identity(2)) + 0.5 * tensor(identity(N_cutoff), sigmaz()))
Hp = g*(tensor(destroy(N_cutoff),sigmap()) + tensor(create(N_cutoff), sigmam()))
H0 + Hp



    


















    
Out[15]:







$${\left( {{\hat{{n}}_{{40}}}\tiny\otimes\normalsize{\hat{{I}}_{{2}}}}+{{0.5}\tiny\times\normalsize{{\hat{{I}}_{{40}}}\tiny\otimes\normalsize{\hat{σ}_{z}}}}+{{0.1}\tiny\times\normalsize{\left( {{\hat{{a^\dagger}}_{{40}}}\tiny\otimes\normalsize{\hat{σ}_{-}}}+{{\hat{{a}}_{{40}}}\tiny\otimes\normalsize{\hat{σ}_{+}}} \right)}} \right)}$$

Definite Photon State





In [16]:



    


n = 3
n_p = tensor(basis(N_cutoff,n), basis(2,0))
np1_m = tensor(basis(N_cutoff,n+1), basis(2,1))
n_p



    


















    
Out[16]:







$${|{3}_{\tiny N\normalsize 40},{0}_{\tiny N\normalsize 2}\rangle}$$
















In [17]:



    


res = sesolve(H0 + Hp, n_p, ts)



    


















    






Generating cython code...
Compiling cython code...










    











  
  
  






#cython: boundscheck=False, wraparound=False, initializedcheck=False, cdivision=True
cimport numpy as np
import numpy as np
from libc.math cimport sin, cos, exp, sinh, cosh, tan, tanh, log
from scipy.linalg.cython_blas cimport zgemv, zgemm

cdef int iONE=1, iZERO=0
cdef double complex zONE=1, zZERO=0, zNHALF=-0.5

# Declaration of global variables for
# - intermediate results
# - predefined constants
cdef double complex [:, :] C = np.empty((80, 1), dtype=np.complex)
cdef double complex [:] A # initialized in setup
cdef int [:] A_pointers # initialized in setup
cdef int [:] A_indices # initialized in setup

# The generated cython code
cdef generated_function(
    # input arguments
    double E,
    double complex [:, :] B,
    # output argument
    double complex [:, :] D
    ):
    # indices and intermediate values for various matrix operations
    cdef int i, j, k, ii, jj, kk
    cdef double complex c
    # intermediate results
    global C
    # predefined constants
    global A, A_pointers, A_indices
    # evaluating the function
    # C = A.B
    ii = A_pointers.shape[0]-1
    for i in range(ii):
        jj = A_pointers[i]
        kk = A_pointers[i+1]
        c = 0
        for j in range(jj,kk):
            k = A_indices[j]
            c += A[j]*B[k,0]
        C[i,0] = c
    
    # D = -1j*C
    ii = C.shape[0]
    jj = C.shape[1]
    for i in range(ii):
        for j in range(jj):
            D[i,j] = -1j*C[i,j]

# Function to deliver the constants from python to cython
cpdef setup_generated_function(
    double complex [:] _A,
    int [:] _A_pointers,
    int [:] _A_indices
    ):
    global A, A_pointers, A_indices
    A = _A
    A_pointers = _A_pointers
    A_indices = _A_indices

ctypedef void rhs_t(double t, double *y, double *ydot)
cdef extern from "cvode_simple_interface.h":
    struct cvsi_instance:
        pass
    cvsi_instance *cvsi_setup(rhs_t *rhs,
                              double *y0, int neq,
                              double t0,
                              long mxsteps, double reltol, double abstol)
    int cvsi_step(cvsi_instance *instance, double t)
    void cvsi_destroy(cvsi_instance *instance)

cdef inline void RHS(double t, double *y, double *ydot):
    generated_function(t,
                       <double complex [:80, :1]> <double complex *> y,
                       <double complex [:80, :1]> <double complex *> ydot)
    # TODO The casts above need to allocate new memview structs each time
    # It is noticeable on 2x2 matrices (10%), not noticeable otherwise

cpdef list pythonsolve(
        np.ndarray[np.double_t, ndim=1] ts,
        np.ndarray[np.complex_t, ndim=2] y0,
        long mxsteps, double rtol, double atol,
        progressbar):
    if ts[0] == 0:
        ts = ts[1:]
        res = [y0]
    else:
        res = []
    cdef np.ndarray[np.complex_t, ndim=2, mode="c"] _y0 = np.copy(y0)
    cdef int neq = np.size(_y0)*2
    cdef cvsi_instance *instance = cvsi_setup(RHS, <double *> &_y0[0,0], neq,
                                              0, mxsteps, rtol, atol)
    for t in ts:
        cvsi_step(instance, t)
        res.append(np.copy(_y0))
        progressbar.step()
    cvsi_destroy(instance)

    progressbar.stop()
    return res














    






Running cython code...
Starting at 10/20 13:53:46.
Finishing at 10/20 13:53:46.
Total time: 0 seconds.
Formatting the output...

















In [18]:



    


ovlps = overlap([n_p, np1_m], res)
plt.plot(ts*g/np.pi, np.abs(ovlps)**2)
plt.legend([r'$|%d,+\rangle$'%n, r'$|%d,-\rangle$'%(n+1)])
plt.title(r'Population-vs-$gt/\pi$');



    


















    
























In [19]:



    


n = 8
n_p = tensor(basis(N_cutoff,n), basis(2,0))
np1_m = tensor(basis(N_cutoff,n+1), basis(2,1))
res = sesolve(H0 + Hp, n_p, ts)
ovlps = overlap([n_p, np1_m], res)
plt.plot(ts*g/np.pi, np.abs(ovlps)**2)
plt.legend([r'$|%d,+\rangle$'%n, r'$|%d,-\rangle$'%(n+1)])
plt.title(r'Population-vs-$gt/\pi$');



    


















    






Generating cython code...
Compiling cython code...










    











  
  
  






#cython: boundscheck=False, wraparound=False, initializedcheck=False, cdivision=True
cimport numpy as np
import numpy as np
from libc.math cimport sin, cos, exp, sinh, cosh, tan, tanh, log
from scipy.linalg.cython_blas cimport zgemv, zgemm

cdef int iONE=1, iZERO=0
cdef double complex zONE=1, zZERO=0, zNHALF=-0.5

# Declaration of global variables for
# - intermediate results
# - predefined constants
cdef double complex [:, :] C = np.empty((80, 1), dtype=np.complex)
cdef double complex [:] A # initialized in setup
cdef int [:] A_pointers # initialized in setup
cdef int [:] A_indices # initialized in setup

# The generated cython code
cdef generated_function(
    # input arguments
    double E,
    double complex [:, :] B,
    # output argument
    double complex [:, :] D
    ):
    # indices and intermediate values for various matrix operations
    cdef int i, j, k, ii, jj, kk
    cdef double complex c
    # intermediate results
    global C
    # predefined constants
    global A, A_pointers, A_indices
    # evaluating the function
    # C = A.B
    ii = A_pointers.shape[0]-1
    for i in range(ii):
        jj = A_pointers[i]
        kk = A_pointers[i+1]
        c = 0
        for j in range(jj,kk):
            k = A_indices[j]
            c += A[j]*B[k,0]
        C[i,0] = c
    
    # D = -1j*C
    ii = C.shape[0]
    jj = C.shape[1]
    for i in range(ii):
        for j in range(jj):
            D[i,j] = -1j*C[i,j]

# Function to deliver the constants from python to cython
cpdef setup_generated_function(
    double complex [:] _A,
    int [:] _A_pointers,
    int [:] _A_indices
    ):
    global A, A_pointers, A_indices
    A = _A
    A_pointers = _A_pointers
    A_indices = _A_indices

ctypedef void rhs_t(double t, double *y, double *ydot)
cdef extern from "cvode_simple_interface.h":
    struct cvsi_instance:
        pass
    cvsi_instance *cvsi_setup(rhs_t *rhs,
                              double *y0, int neq,
                              double t0,
                              long mxsteps, double reltol, double abstol)
    int cvsi_step(cvsi_instance *instance, double t)
    void cvsi_destroy(cvsi_instance *instance)

cdef inline void RHS(double t, double *y, double *ydot):
    generated_function(t,
                       <double complex [:80, :1]> <double complex *> y,
                       <double complex [:80, :1]> <double complex *> ydot)
    # TODO The casts above need to allocate new memview structs each time
    # It is noticeable on 2x2 matrices (10%), not noticeable otherwise

cpdef list pythonsolve(
        np.ndarray[np.double_t, ndim=1] ts,
        np.ndarray[np.complex_t, ndim=2] y0,
        long mxsteps, double rtol, double atol,
        progressbar):
    if ts[0] == 0:
        ts = ts[1:]
        res = [y0]
    else:
        res = []
    cdef np.ndarray[np.complex_t, ndim=2, mode="c"] _y0 = np.copy(y0)
    cdef int neq = np.size(_y0)*2
    cdef cvsi_instance *instance = cvsi_setup(RHS, <double *> &_y0[0,0], neq,
                                              0, mxsteps, rtol, atol)
    for t in ts:
        cvsi_step(instance, t)
        res.append(np.copy(_y0))
        progressbar.step()
    cvsi_destroy(instance)

    progressbar.stop()
    return res














    






Running cython code...
Starting at 10/20 13:53:58.
Finishing at 10/20 13:53:58.
Total time: 0 seconds.
Formatting the output...

Coherent State





In [20]:



    


alpha = 5
coh = tensor(coherent(N_cutoff, alpha), basis(2,0))
coh



    


















    
Out[20]:







$${|{\tiny\alpha\normalsize 5.00}_{\tiny N\normalsize 40},{0}_{\tiny N\normalsize 2}\rangle}$$
















In [21]:



    


ts = 80/g*np.linspace(0,1,4000)
res = sesolve(H0 + Hp, coh, ts)
inversion = expect(tensor(identity(N_cutoff), sigmaz()), res)



    


















    






Generating cython code...
Compiling cython code...










    











  
  
  






#cython: boundscheck=False, wraparound=False, initializedcheck=False, cdivision=True
cimport numpy as np
import numpy as np
from libc.math cimport sin, cos, exp, sinh, cosh, tan, tanh, log
from scipy.linalg.cython_blas cimport zgemv, zgemm

cdef int iONE=1, iZERO=0
cdef double complex zONE=1, zZERO=0, zNHALF=-0.5

# Declaration of global variables for
# - intermediate results
# - predefined constants
cdef double complex [:, :] C = np.empty((80, 1), dtype=np.complex)
cdef double complex [:] A # initialized in setup
cdef int [:] A_pointers # initialized in setup
cdef int [:] A_indices # initialized in setup

# The generated cython code
cdef generated_function(
    # input arguments
    double E,
    double complex [:, :] B,
    # output argument
    double complex [:, :] D
    ):
    # indices and intermediate values for various matrix operations
    cdef int i, j, k, ii, jj, kk
    cdef double complex c
    # intermediate results
    global C
    # predefined constants
    global A, A_pointers, A_indices
    # evaluating the function
    # C = A.B
    ii = A_pointers.shape[0]-1
    for i in range(ii):
        jj = A_pointers[i]
        kk = A_pointers[i+1]
        c = 0
        for j in range(jj,kk):
            k = A_indices[j]
            c += A[j]*B[k,0]
        C[i,0] = c
    
    # D = -1j*C
    ii = C.shape[0]
    jj = C.shape[1]
    for i in range(ii):
        for j in range(jj):
            D[i,j] = -1j*C[i,j]

# Function to deliver the constants from python to cython
cpdef setup_generated_function(
    double complex [:] _A,
    int [:] _A_pointers,
    int [:] _A_indices
    ):
    global A, A_pointers, A_indices
    A = _A
    A_pointers = _A_pointers
    A_indices = _A_indices

ctypedef void rhs_t(double t, double *y, double *ydot)
cdef extern from "cvode_simple_interface.h":
    struct cvsi_instance:
        pass
    cvsi_instance *cvsi_setup(rhs_t *rhs,
                              double *y0, int neq,
                              double t0,
                              long mxsteps, double reltol, double abstol)
    int cvsi_step(cvsi_instance *instance, double t)
    void cvsi_destroy(cvsi_instance *instance)

cdef inline void RHS(double t, double *y, double *ydot):
    generated_function(t,
                       <double complex [:80, :1]> <double complex *> y,
                       <double complex [:80, :1]> <double complex *> ydot)
    # TODO The casts above need to allocate new memview structs each time
    # It is noticeable on 2x2 matrices (10%), not noticeable otherwise

cpdef list pythonsolve(
        np.ndarray[np.double_t, ndim=1] ts,
        np.ndarray[np.complex_t, ndim=2] y0,
        long mxsteps, double rtol, double atol,
        progressbar):
    if ts[0] == 0:
        ts = ts[1:]
        res = [y0]
    else:
        res = []
    cdef np.ndarray[np.complex_t, ndim=2, mode="c"] _y0 = np.copy(y0)
    cdef int neq = np.size(_y0)*2
    cdef cvsi_instance *instance = cvsi_setup(RHS, <double *> &_y0[0,0], neq,
                                              0, mxsteps, rtol, atol)
    for t in ts:
        cvsi_step(instance, t)
        res.append(np.copy(_y0))
        progressbar.step()
    cvsi_destroy(instance)

    progressbar.stop()
    return res














    






Running cython code...
Starting at 10/20 13:54:08.
Finishing at 10/20 13:54:13.
Total time: 5 seconds.
Formatting the output...

















In [22]:



    


plt.plot(ts*g, inversion)
plt.title(r'$\langle \hat{\sigma}_z \rangle$-vs-$gt$');

Content source: Krastanov/cutiepy

Similar notebooks:

notebook.community | gallery | about