我试图调整 qutip Wigner 函数，让它处理两种模式状态，特别是迭代方法。
但是我的输出给出的数组大小太大，我不确定为什么？那就是当我尝试使用它计算 Wigner 对数负性时，积分以数组而不是奇异值的形式出现。
它的用途的代码和描述如下：
`import numpy as np
from scipy import (zeros, array, arange, exp, real, conj, pi,
copy, sqrt, meshgrid, size, polyval, fliplr, conjugate,
cos, sin)
import scipy.sparse as sp
import scipy.fftpack as ft
import scipy.linalg as la
from scipy.special import genlaguerre
from scipy.special import binom
from scipy.special import sph_harm
from qutip.qobj import Qobj, isket, isoper
from qutip.states import ket2dm
from qutip.parallel import parfor
from qutip.utilities import clebsch
from scipy.special import factorial
from qutip.cy.sparse_utils import _csr_get_diag
from qutip import *
def wigner2(psi, xvec1, yvec1, xvec2, yvec2, method='iterative', g=np.sqrt(2)):
"""Wigner function for a state vector or density matrix at points
`xvec1 + i * yvec1` `xvec2 + i * yvec2`
参数
state : qobj 状态向量或密度矩阵。
xvec1 : array_like x 坐标，用于计算 Wigner 函数。
yvec1 : array_like 计算 Wigner 函数的 y 坐标。
xvec2 : array_like x 坐标，用于计算 Wigner 函数。
yvec2 : array_like 计算 Wigner 函数的 y 坐标。
g: float 缩放因子a = 0.5 * g * (x + iy)，默认g = sqrt(2)。
方法：字符串 {'iterative'} 选择方法 'iterative'，其中 'iterative' 使用迭代方法来评估密度矩阵的 Wigner 函数 :math: |m><n|。'iterative' 方法是默认的，通常推荐使用，但'laguerre' 方法对于非常稀疏的密度矩阵（例如，大希尔伯特空间中的 Fock 状态的叠加）更有效。'fft' 方法是处理具有大量激发 (>~50) 的密度矩阵的首选方法。
退货
W : 表示在指定范围 [xvec1,yvec1] 和 [xvec2,yvec2] """ 上计算的 Wigner 函数的数组值
if not (psi.type == 'ket' or psi.type == 'oper' or psi.type == 'bra'):
raise TypeError('Input state is not a valid operator.')
if psi.type == 'ket' or psi.type == 'bra':
rho = ket2dm(psi) #always use density matrix
else:
rho = psi
if method == 'iterative':
return _wigner2_iterative(rho, xvec1, yvec1, xvec2, yvec2, g)
else:
raise TypeError(
"method must be 'iterative'")
def _wigner2_iterative(rho, xvec1, yvec1, xvec2, yvec2, g=np.sqrt(2)):
"""使用迭代方法评估 Fock 状态的 wigner 函数 :math: |mp><nq|
Wigner 函数计算为 :math:W = \sum_{mpnq} \\rho_{mpnq} W_{mpnq}其中 :math:W_{mpnq}是密度矩阵 :math: 的 Wigner 函数|mp><nq|。在此实现中，对于每个第 m*p 行，Wlist 包含 Wigner 函数 Wlist = [0, ..., W_mpmp, ..., W_mpnq]。只要计算出一个 W_mpnq Wigner 函数，相应的贡献就会添加到总 Wigner 函数中，加权通过密度矩阵中的相应元素 :math: rho_{mpnq}."""
M1 = np.prod(ptrace(rho, 0).shape[0])
M2 = np.prod(ptrace(rho, 1).shape[0])
M = np.prod(rho.shape[0])
X1, Y1, X2, Y2 = np.meshgrid(xvec1, yvec1, xvec2, yvec2)
A1 = 0.5 * g * (X1 + 1.0j * Y1 + 0 * X2 + 0 * Y2)
A2 = 0.5 * g * (0 * X1 + 0 * Y1 + X2 + 1.0j * Y2)
Wlist1 = array([zeros(np.shape(A1), dtype=complex) for k in range(M)])
Wlist2 = array([zeros(np.shape(A2), dtype=complex) for k in range(M)])
W = real(rho[0, 0]) * real(Wlist1[0] * Wlist2[0])
for m in range(0,M1):
if m==0:
Wlist1[0] = exp(-2.0 * abs(A1) ** 2) / (pi)
else:
Wlist1[m] = ((2.0 * A1 * Wlist1[m - 1]) / sqrt(m))
for n in range(0, M2):
if n==0:
Wlist2[0] = exp(-2.0 * abs(A2) ** 2) / (pi)
else:
Wlist2[n] = ((2.0 * A2 * Wlist2[n - 1]) / sqrt(n))
if m != 0 and n != 0:
W += 2 * real(rho[0, m * M2 + n] * Wlist1[m] * Wlist2[n])
for p in range(0, M1):
temp1 = copy(Wlist1[m])
temp2 = copy(Wlist2[n])
if p==0:
Wlist1[p] = exp(-2.0 * abs(A1) ** 2) / (pi)
else:
Wlist1[p] = ((2.0 * conj(A1) * temp1 -sqrt(p) * Wlist1[p-1]) / sqrt(p))
for q in range(0, M2):
if q==0:
Wlist2[q] = exp(-2.0 * abs(A2) ** 2) / (pi)
else:
Wlist2[q] = ((2.0 * conj(A2) * temp2 - sqrt(q) * Wlist1[q - 1]) / sqrt(q))
W += 2 * real(rho[p * M2 + q, p * M2 + q] * Wlist1[p] * Wlist2[q])
if p != 0 and q !=0:
for k in range(p + 1, M1):
temp3 = (2 * A1 * Wlist1[k-1] - sqrt(k) * temp1) / sqrt(k)
temp1 = copy(Wlist1[k])
Wlist1[k] = temp3
for l in range(q +1, M2):
temp4 = (2 * A2 * Wlist2[l-1] - sqrt(l) * temp2) / sqrt(l)
temp2 = copy(Wlist2[l])
Wlist2[l] = temp4
W += 2 * real(rho[p * M2 + q, k *M2 +l] * Wlist1[k] * Wlist2[l])
return 0.5 * W * g **2'