Robert Johansson, http://jrjohansson.github.com, robert@riken.jp
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from wavefunction import *
from wavefunction.wavefunction2d import *
Here follow Orlando et al., PRB 60, 15398 (1999) and numerically calculate the wavefunctions and energy-levels of the Flux qubit device using a two-dimensional potential of the superconducting Flux qubit circuit. The two generalized coordinates are the gauge invariant phases across the Josephson junctions in the device.
The Hamiltonian in the original phase coordinates is
$H_t = \frac{1}{2} P^T M^{-1} P + E_J(2 + \alpha - \cos\varphi_1 - \cos\varphi_2 - \alpha\cos(2\pi f + \varphi_1 - \varphi_2))$
where $P = (-i\hbar\partial_{\varphi_1}, -i\hbar\partial_{\varphi_2})^T$ is the generalized momentum vector and the mass matrix is
$M = (\Phi_0/2\pi)^2 C \begin{pmatrix} 1 + \alpha + \gamma & -\alpha \\ -\alpha & 1 + \alpha + \gamma\end{pmatrix}$
and
$M^{-1} = \frac{1}{(\Phi_0/2\pi)^2 C} \frac{1}{(1 + \gamma)(1 + \gamma + 2\alpha)} \begin{pmatrix} 1 + \alpha + \gamma & \alpha \\ \alpha & 1 + \alpha + \gamma\end{pmatrix}$
This gives
$H_t = - \frac{4E_C }{(1 + \gamma)(1 + \gamma + 2\alpha)} \left[ (1 + \gamma + \alpha)\left(\frac{\partial}{\partial\varphi_1}\right)^2 + 2\alpha\frac{\partial}{\partial\varphi_1}\frac{\partial}{\partial\varphi_2} + (1 + \gamma + \alpha)\left(\frac{\partial}{\partial\varphi_2}\right)^2 \right]
For the numerical calculations we use units of $E_J$
$H_t = - \frac{4}{(1 + \gamma)(1 + \gamma + 2\alpha)}\frac{E_C}{E_J} \left[ (1 + \gamma + \alpha)\left(\frac{\partial}{\partial\varphi_1}\right)^2 + 2\alpha\frac{\partial}{\partial\varphi_1}\frac{\partial}{\partial\varphi_2} + (1 + \gamma + \alpha)\left(\frac{\partial}{\partial\varphi_2}\right)^2 \right]
args = {'Ec': 1/80.0, 'Ej': 1.0, 'alpha': 0.8, 'gamma': 0.0, 'f': 0.50}
x1 = np.linspace(-2*np.pi, 2*np.pi, 100)
x2 = np.linspace(-2*np.pi, 2*np.pi, 100)
X1, X2 = np.meshgrid(x1, x2)
def U_flux_qubit(x1, x2, args):
globals().update(args)
return 2 + alpha - np.cos(x1) - np.cos(x2) - alpha * np.cos(2 * np.pi * f + x1 - x2)
Z = np.real(U_flux_qubit(X1, X2, args))
fig, ax = plt.subplots()
p = ax.pcolor(X1/(2*np.pi), X2/(2*np.pi), Z, cmap=cm.RdBu, vmin=Z.min(), vmax=Z.max())
ax.set_xlabel(r'$x_1$', fontsize=16)
ax.set_ylabel(r'$x_2$', fontsize=16)
cb = fig.colorbar(p, ax=ax)
$H_t = - \frac{4}{(1 + \gamma)(1 + \gamma + 2\alpha)}\frac{E_C}{E_J} \left[ (1 + \gamma + \alpha)\left(\frac{\partial}{\partial\varphi_1}\right)^2 + 2\alpha\frac{\partial}{\partial\varphi_1}\frac{\partial}{\partial\varphi_2} + (1 + \gamma + \alpha)\left(\frac{\partial}{\partial\varphi_2}\right)^2 \right]
$\displaystyle K_{n_1, n_2}^{m_1, m_2} = \delta_{n_1,m_1} \delta_{n_2,m_2} \left(k_{11}\left(\frac{2\pi m_1}{T_{x_1}}\right)^2 + k_{12}\frac{(2\pi)^2 m_1m_2}{T_{x_1}T_{x_2}} + k_{22}\left(\frac{2\pi m_2}{T_{x_2}}\right)^2\right)$
# pick a truncation of the fourier series: n1 = [-L1, ..., L1], n2 = [-L2, ..., L2]
L1 = 5
L2 = 5
# pick periods for the coordinates
Tx1 = Tx2 = 2 * np.pi
#
k11 = k22 = 4 * (Ec / Ej) * (1 + alpha + gamma) / ((1 + gamma) * (1 + 2 * alpha + gamma))
k12 = 4 * (Ec / Ej) * 2 * alpha / ((1 + gamma) * (1 + 2 * alpha + gamma))
K = assemble_K(L1, L2, k11, k12, k22, Tx1, Tx2)
The flux qubit potential we consider here is
$\displaystyle U(\varphi_1, \varphi_2)/E_J = 2 + \alpha - \cos(\varphi_1) - \cos(\varphi_2) - \alpha \cos(2\pi f + \varphi_1 - \varphi_2)$
To obtain the Fourier series expansion of $U(\varphi_1, \varphi_2)$, we first write the $\cos$ expressions as exponential functions
$\displaystyle U(\varphi_1, \varphi_2)/E_J = (2 + \alpha)
$\displaystyle U(\varphi_1, \varphi_2)/E_J = \sum_{n_1}\sum_{n_2} u_{n_1n_2} e^{in_1\varphi_1}e^{in_2\varphi_2}$
we can identity
$\displaystyle u_{n_1n_2} = (2 + \alpha) \delta_{0,n_1}\delta_{0,n_2}
def assemble_u_flux_qubit(L1, L2, args, sparse=False):
globals().update(args)
L1n = 2 * L1 + 1
L2n = 2 * L2 + 1
u = np.zeros((L1n, L2n), dtype=np.complex)
for n1 in np.arange(-L1, L1+1):
N1 = n1 + L1
for n2 in np.arange(-L2, L2+1):
N2 = n2 + L2
u[N1, N2] = (2 + alpha) * delta(0, n1) * delta(0, n2) + \
- 0.5 * (delta(1, n1) + delta(-1, n1)) * delta(0, n2) + \
- 0.5 * (delta(1, n2) + delta(-1, n2)) * delta(0, n1)+ \
- 0.5 * alpha * np.exp(+2j * np.pi * f) * delta(+1, n1) * delta(-1, n2) + \
- 0.5 * alpha * np.exp(-2j * np.pi * f) * delta(-1, n1) * delta(+1, n2)
return u
u = assemble_u_flux_qubit(L1, L2, args)
$\displaystyle U(x_p, x_m) = \sum_{n_1}\sum_{n_2} u_{n_1n_2} e^{in_1x_p}e^{in_2x_m}$
Z = np.real(evalute_fourier_series(X1, X2, L1, L2, u))
fig, ax = plt.subplots()
p = ax.pcolor(X1/(2*np.pi), X2/(2*np.pi), Z, cmap=cm.RdBu, vmin=Z.min(), vmax=Z.max())
cb = fig.colorbar(p, ax=ax)
Given this Fourier series we can calcuate the potential contribution to the eigenvalue problem as
$\displaystyle V_{n_1, n_2}^{m_1, m_2} = u_{n_1-m_1,n_2-m_2}$
V = assemble_V(L1, L2, u)
and specifically for the flux qubit potential
$\displaystyle V_{n_1, n_2}^{m_1, m_2} = (2 + \alpha) \delta_{0,n_1-m_1}\delta_{0,n_2-m_2}
def assemble_V_flux_qubit(L1, L2, args, sparse=False):
globals().update(args)
L1n = 2 * L1 + 1
L2n = 2 * L2 + 1
V = np.zeros((L1n*L1n, L2n*L2n), dtype=np.complex)
for n1 in np.arange(-L1, L1+1):
for n2 in np.arange(-L1, L1+1):
N = index_m2v(L1, n1, n2)
for m1 in np.arange(-L2, L2+1):
for m2 in np.arange(-L2, L2+1):
M = index_m2v(L2, m1, m2)
V[N,M] = (2 + alpha) * delta(m1, n1) * delta(m2, n2) + \
- 0.5 * (delta(m1 + 1, n1) + delta(m1 - 1, n1)) * delta(m2, n2) + \
- 0.5 * (delta(m2 + 1, n2) + delta(m2 - 1, n2)) * delta(m1, n1) + \
- 0.5 * alpha * np.exp(+2j * np.pi * f) * delta(m1 + 1, n1) * delta(m2 - 1, n2) + \
- 0.5 * alpha * np.exp(-2j * np.pi * f) * delta(m1 - 1, n1) * delta(m2 + 1, n2)
return V
V_fq = assemble_V_flux_qubit(L1, L2, args)
abs(V - V_fq).max()
0.0
We now want to solve the eigenvalue problem
$H \Psi = E \Psi$
where $H$ is the matrix representation of the Hamiltonian assembled from $K$ and $V$ above.
H = K + V
vals, vecs = solve_eigenproblem(H)
fig, ax = plt.subplots()
U_x = np.real(U_flux_qubit(x1, -x1, args))
for val in vals:
ax.plot(x1, np.real(val) * np.ones_like(x1), 'k')
ax.plot(x1, U_x, label="potential", lw='2')
ax.axis('tight')
ax.set_ylim(min(vals[0], U_x.min()), U_x.max())
ax.set_ylabel(r'$U(x_1, x_2=x_1)$', fontsize=16)
ax.set_xlabel(r'$x_1$', fontsize=16);
Nstates = 4
fig, axes = plt.subplots(Nstates, 3, figsize=(16, 5 * Nstates))
for n in range(Nstates):
psi = convert_v2m(L1, L2, vecs[n])
Z = np.real(evalute_fourier_series(X1, X2, L1, L2, u))
PSI = evalute_fourier_series(X1, X2, L1, L2, psi)
p = axes[n, 0].pcolor(X1/(2*np.pi), X2/(2*np.pi), np.real(PSI),
cmap=cm.RdBu, vmin=-abs(PSI.real).max(), vmax=abs(PSI.real).max())
c = axes[n, 0].contour(X1/(2*np.pi), X2/(2*np.pi), Z, cmap=cm.RdBu, vmin=Z.min(), vmax=Z.max())
cb = fig.colorbar(p, ax=axes[n, 0])
p = axes[n, 1].pcolor(X1/(2*np.pi), X2/(2*np.pi), np.imag(PSI),
cmap=cm.RdBu, vmin=-abs(PSI.imag).max(), vmax=abs(PSI.imag).max())
c = axes[n, 1].contour(X1/(2*np.pi), X2/(2*np.pi), Z, cmap=cm.RdBu, vmin=Z.min(), vmax=Z.max())
cb = fig.colorbar(p, ax=axes[n, 1])
p = axes[n, 2].pcolor(X1/(2*np.pi), X2/(2*np.pi), np.abs(PSI), cmap=cm.Blues, vmin=0, vmax=abs(PSI).max())
c = axes[n, 2].contour(X1/(2*np.pi), X2/(2*np.pi), Z, cmap=cm.RdBu, vmin=Z.min(), vmax=Z.max())
cb = fig.colorbar(p, ax=axes[n, 2])
axes[n, 0].set_ylabel("%d state" % n);
axes[0, 0].set_title(r"$\mathrm{re}\;\Psi(\varphi_p, \varphi_m)$", fontsize=16);
axes[0, 1].set_title(r"$\mathrm{im}\;\Psi(\varphi_p, \varphi_m)$", fontsize=16);
axes[0, 2].set_title(r"$|\Psi(\varphi_p, \varphi_m)|^2$", fontsize=16);
f_vec = np.linspace(0.45, 0.55, 50)
e_vals = np.zeros((len(vals), len(f_vec)))
K = assemble_K(L1, L2, k11, k12, k22, Tx1, Tx2)
for f_idx, f in enumerate(f_vec):
args['f'] = f
#u = assemble_u_flux_qubit(L1, L2, args)
#V = assemble_V(L1, L2, u)
V = assemble_V_flux_qubit(L1, L2, args)
H = K + V
vals, vecs = solve_eigenproblem(H)
e_vals[:, f_idx] = np.real(vals)
fig, ax = plt.subplots(1, 1, figsize=(10,6))
for n in range(len(vals)):
ax.plot(f_vec, e_vals[n, :])
ax.axis('tight')
ax.set_ylim(e_vals[0, :].min(), e_vals[20, :].max())
ax.set_ylabel("Energy levels")
ax.set_xlabel(r'$f$');
The Hamiltonian in the rotated coordinates, where $\phi_p = (\phi_1 + \phi_2)/2$ and $\phi_m = (\phi_1 - \phi_2)/2$, is
$H_t = \frac{1}{2M_p}\left(-i\hbar\frac{\partial}{\partial\phi_p}\right)^2 + \frac{1}{2M_m}\left(-i\hbar\frac{\partial}{\partial\phi_m}\right)^2 + E_J(2 + \alpha - 2 \cos\phi_p\cos\phi_m - \alpha\cos(2\pi f + 2 \phi_m))$
where
$M_p = \left(\frac{\Phi_0}{2\pi}\right)^2 2 C(1 + \gamma),$
$M_m = \left(\frac{\Phi_0}{2\pi}\right)^2 2 C(1 + 2\alpha + \gamma).$
For the numerical calculations we use units of $E_J$
$H_t = -\frac{2}{1 + \gamma}\frac{E_C}{E_J}\left(\frac{\partial}{\partial\phi_p}\right)^2
args = {'Ec': 1/80.0, 'Ej': 1.0, 'alpha': 0.8, 'gamma': 0.0, 'f': 0.50}
x_p = np.linspace(-2*np.pi, 2*np.pi, 100)
x_m = np.linspace(-2*np.pi, 2*np.pi, 100)
X1, X2 = np.meshgrid(x_p, x_m)
def U_flux_qubit(x_p, x_m, args):
globals().update(args)
return 2 + alpha - 2 * np.cos(x_p) * np.cos(x_m) - alpha * np.cos(2 * np.pi * f + 2 * x_m)
Z = np.real(U_flux_qubit(X1, X2, args))
fig, ax = plt.subplots()
p = ax.pcolor(X1/(2*np.pi), X2/(2*np.pi), Z, cmap=cm.RdBu, vmin=Z.min(), vmax=Z.max())
ax.set_xlabel(r'$x_p$', fontsize=16)
ax.set_ylabel(r'$x_m$', fontsize=16)
cb = fig.colorbar(p, ax=ax)
$\displaystyle K_{n_1, n_2}^{m_1, m_2} = \delta_{n_1,m_1} \delta_{n_2,m_2} \left(\hbar^2 k_{11}\left(\frac{2\pi m_1}{T_{x_1}}\right)^2 + \hbar^2 k_{12}\frac{(2\pi)^2 m_1m_2}{T_{x_1}T_{x_2}} + \hbar^2 k_{22}\left(\frac{2\pi m_2}{T_{x_2}}\right)^2\right)$
# pick a truncation of the fourier series: n1 = [-L1, ..., L1], n2 = [-L2, ..., L2]
L1 = 10
L2 = 10
# pick periods for the coordinates
Tx1 = Tx2 = 2 * np.pi
#
k11, k12, k22 = 2 / (1 + gamma) * Ec / Ej, 0.0, 2 / (1 + 2 * alpha + gamma) * Ec / Ej
K = assemble_K(L1, L2, k11, k12, k22, Tx1, Tx2)
The potential contribution to the matrix representation of the Hamiltonian can be calculated analytically or numerically. Here we first focus on the analytical calculation for a flux qubit potential (period with period $2\pi$), and then discuss how it could be evaluated completely numerically.
The flux qubit potential we consider here is
$\displaystyle U(x_p, x_m) = 2 + \alpha - 2 \cos\phi_p\cos\phi_m - \alpha\cos(2\pi f + 2 \phi_m)$
To obtain the Fourier series expansion of $U(x_p, x_m)$, we first write the $\cos$ expressions as exponential functions
$\displaystyle U(x_p, x_m) = 2 + \alpha
Now, by comparing with the Fourier series
$\displaystyle U(x_p, x_m) = \sum_{n_1}\sum_{n_2} u_{n_1n_2} e^{in_1x_p}e^{in_2x_m}$
we can identity
$\displaystyle u_{n_1n_2} = (2 + \alpha) \delta_{0,n_1}\delta_{0,n_2}
\frac{1}{2} ( \delta_{1, n_1}\delta_{1, n_2} + \delta_{1, n_1}\delta_{-1,n_2} + \delta_{-1, n_1}\delta_{1,n_2} + \delta_{-1, n_1}\delta_{-1,n_2} )
\frac{\alpha}{2} (e^{i2\pi f} \delta_{0,n_1}\delta_{+2,n_2} + e^{-i2\pi f}\delta_{0,n_1}\delta_{-2,n_2}) $
def assemble_u_flux_qubit(L1, L2, args, sparse=False):
globals().update(args)
L1n = 2 * L1 + 1
L2n = 2 * L2 + 1
u = np.zeros((L1n, L2n), dtype=np.complex)
for n1 in np.arange(-L1, L1+1):
N1 = n1 + L1
for n2 in np.arange(-L2, L2+1):
N2 = n2 + L2
u[N1, N2] = (2 + alpha) * delta(0, n1) * delta(0, n2) + \
- 0.5 * (delta(+1, n1) * delta(+1, n2) + delta(+1, n1) * delta(-1, n2) +
delta(-1, n1) * delta(+1, n2) + delta(-1, n1) * delta(-1, n2)) + \
- 0.5 * alpha * np.exp(+2j * np.pi * f) * delta(0, n1) * delta(+2, n2) + \
- 0.5 * alpha * np.exp(-2j * np.pi * f) * delta(0, n1) * delta(-2, n2)
return u
def assemble_u_flux_qubit_direct(L1, L2, args, sparse=False):
globals().update(args)
L1n = 2 * L1 + 1
L2n = 2 * L2 + 1
u = np.zeros((L1n, L2n), dtype=np.complex)
u[L1, L2] = 2 + alpha
u[+1+L1, +1+L2] = - 0.5
u[-1+L1, +1+L2] = - 0.5
u[+1+L1, -1+L2] = - 0.5
u[-1+L1, -1+L2] = - 0.5
u[L1, +2+L2] = -0.5 * alpha * np.exp(+2j * np.pi * f)
u[L1, -2+L2] = -0.5 * alpha * np.exp(-2j * np.pi * f)
return u
u = assemble_u_flux_qubit(L1, L2, args)
u_direct = assemble_u_flux_qubit_direct(L1, L2, args)
abs(u - u_direct).max() == 0
True
$\displaystyle U(x_p, x_m) = \sum_{n_1}\sum_{n_2} u_{n_1n_2} e^{in_1x_p}e^{in_2x_m}$
Z = np.real(evalute_fourier_series(X1, X2, L1, L2, u))
fig, ax = plt.subplots()
p = ax.pcolor(X1/(2*np.pi), X2/(2*np.pi), Z, cmap=cm.RdBu, vmin=Z.min(), vmax=Z.max())
ax.set_xlabel(r'$x_p$', fontsize=16)
ax.set_ylabel(r'$x_m$', fontsize=16)
cb = fig.colorbar(p, ax=ax)
Given this Fourier series we can calcuate the potential contribution to the eigenvalue problem as
$\displaystyle V_{n_1, n_2}^{m_1, m_2} = u_{n_1-m_1,n_2-m_2}$
and specifically for the flux qubit potential
$\displaystyle V_{n_1, n_2}^{m_1, m_2} = (2 + \alpha) \delta_{0,n_1-m_1}\delta_{0,n_2-m_2}
\frac{1}{2} ( \delta_{1, n_1-m_1}\delta_{1, n_2-m_2} + \delta_{1, n_1-m_1}\delta_{-1,n_2-m_2} + \delta_{-1, n_1-m_1}\delta_{1,n_2-m_2} + \delta_{-1, n_1-m_1}\delta_{-1,n_2-m_2} )
\frac{\alpha}{2} (e^{i2\pi f} \delta_{0,n_1-m_1}\delta_{+2,n_2-m_2} + e^{-i2\pi f}\delta_{0,n_1-m_1}\delta_{-2,n_2-m_2}) $
def assemble_V_flux_qubit(L1, L2, args, sparse=False):
globals().update(args)
L1n = 2 * L1 + 1
L2n = 2 * L2 + 1
V = np.zeros((L1n*L1n, L2n*L2n), dtype=np.complex)
for n1 in np.arange(-L1, L1+1):
for n2 in np.arange(-L1, L1+1):
N = index_m2v(L1, n1, n2)
for m1 in np.arange(-L2, L2+1):
for m2 in np.arange(-L2, L2+1):
M = index_m2v(L2, m1, m2)
V[N,M] = (2 + alpha) * delta(m1, n1) * delta(m2, n2) + \
- 0.5 * (delta(m1 + 1, n1) * delta(m2 + 1, n2) +
delta(m1 + 1, n1) * delta(m2 - 1, n2) +
delta(m1 - 1, n1) * delta(m2 + 1, n2) +
delta(m1 - 1, n1) * delta(m2 - 1, n2)) + \
- 0.5 * alpha * np.exp(+2j * np.pi * f) * delta(m1, n1) * delta(m2 + 2, n2) + \
- 0.5 * alpha * np.exp(-2j * np.pi * f) * delta(m1, n1) * delta(m2 - 2, n2)
return V
V_full = assemble_V(L1, L2, u)
V = assemble_V_flux_qubit(L1, L2, args)
abs(V-V_full).max()
0.0
We now want to solve the eigenvalue problem
$H \Psi = E \Psi$
where $H$ is the matrix representation of the Hamiltonian assembled from $K$ and $V$ above.
H = K + V
vals, vecs = solve_eigenproblem(H)
fig, ax = plt.subplots()
U_x = np.real(U_flux_qubit(0, x_m, args))
for val in vals:
ax.plot(x_m, np.real(val) * np.ones_like(x_m), 'k')
ax.plot(x_m, U_x, label="potential", lw='2')
ax.axis('tight')
ax.set_ylim(min(vals[0], U_x.min()), U_x.max())
ax.set_ylabel(r'$U(x_1, x_2=x_1)$', fontsize=16)
ax.set_xlabel(r'$x_1$', fontsize=16);
Nstates = 6
fig, axes = plt.subplots(Nstates, 3, figsize=(16, 5 * Nstates))
for n in range(Nstates):
psi = convert_v2m(L1, L2, vecs[n])
Z = np.real(evalute_fourier_series(X1, X2, L1, L2, u))
PSI = evalute_fourier_series(X1, X2, L1, L2, psi)
p = axes[n, 0].pcolor(X1/(2*np.pi), X2/(2*np.pi), np.real(PSI),
cmap=cm.RdBu, vmin=-abs(PSI.real).max(), vmax=abs(PSI.real).max())
c = axes[n, 0].contour(X1/(2*np.pi), X2/(2*np.pi), Z, cmap=cm.RdBu, vmin=Z.min(), vmax=Z.max())
cb = fig.colorbar(p, ax=axes[n, 0])
p = axes[n, 1].pcolor(X1/(2*np.pi), X2/(2*np.pi), np.imag(PSI),
cmap=cm.RdBu, vmin=-abs(PSI.imag).max(), vmax=abs(PSI.imag).max())
c = axes[n, 1].contour(X1/(2*np.pi), X2/(2*np.pi), Z, cmap=cm.RdBu, vmin=Z.min(), vmax=Z.max())
cb = fig.colorbar(p, ax=axes[n, 1])
p = axes[n, 2].pcolor(X1/(2*np.pi), X2/(2*np.pi), np.abs(PSI), cmap=cm.Blues, vmin=0, vmax=abs(PSI).max())
c = axes[n, 2].contour(X1/(2*np.pi), X2/(2*np.pi), Z, cmap=cm.RdBu, vmin=Z.min(), vmax=Z.max())
cb = fig.colorbar(p, ax=axes[n, 2])
axes[n, 0].set_ylabel("%d state" % n);
axes[0, 0].set_title(r"$\mathrm{re}\;\Psi(\varphi_p, \varphi_m)$", fontsize=16);
axes[0, 1].set_title(r"$\mathrm{im}\;\Psi(\varphi_p, \varphi_m)$", fontsize=16);
axes[0, 2].set_title(r"$|\Psi(\varphi_p, \varphi_m)|^2$", fontsize=16);
f_vec = np.linspace(0.4, 0.6, 50)
f_vec = np.linspace(0.47, 0.5, 25)
f_vec = np.linspace(0.45, 0.55, 50)
e_vals = np.zeros((len(vals), len(f_vec)))
for f_idx, f in enumerate(f_vec):
args['f'] = f
K = assemble_K(L1, L2, k11, k12, k22, Tx1, Tx2)
V = assemble_V_flux_qubit(L1, L2, args)
H = K + V
vals, vecs = solve_eigenproblem(H)
e_vals[:, f_idx] = np.real(vals)
fig, ax = plt.subplots(1, 1, figsize=(10,6))
for n in range(len(vals)):
ax.plot(f_vec, e_vals[n, :])
ax.axis('tight')
ax.set_ylim(e_vals[0, :].min(), e_vals[20, :].max())
ax.set_ylabel("Energy levels")
ax.set_xlabel(r'$f$');
%reload_ext version_information
%version_information numpy, scipy, matplotlib, wavefunction
Software | Version |
---|---|
Python | 2.7.4 (default, Apr 19 2013, 18:28:01) [GCC 4.7.3] |
IPython | 2.0.0-dev |
OS | posix [linux2] |
numpy | 1.8.0.dev-895866d |
scipy | 0.13.0.dev-7c6d92e |
matplotlib | 1.3.x |
wavefunction | 1.0.0 |
Mon Sep 02 15:30:51 2013 JST |