6

次のようなフォンノイマン方程式があります: dr/dt = - i [H, r]ここで、r と H は複素数の正方行列であり、Python スクリプトを使用して r(t) を見つける必要があります。

そのような方程式を統合するための標準的な手段はありますか?

シュレディンガー方程式 dy/dt = - i H yのように、ベクトルを初期値として別の方程式を解いていたとき、scipy.integrate.ode 関数 ('zvode') を使用しましたが、フォン ノイマンに対して同じ関数を使用しようとしました。方程式は私に次のエラーを与えます:

**scipy/integrate/_ode.py:869: UserWarning: zvode: Illegal input detected. (See printed message.)
ZVODE--  ZWORK length needed, LENZW (=I1), exceeds LZW (=I2)
self.messages.get(istate, 'Unexpected istate=%s' % istate))
  In above message,  I1 =        72   I2 =        24**

コードは次のとおりです。

def integrate(r, t0, t1, dt):
  e = linspace(t0, t1, (t1 - t0) / dt + 10)
  g = linspace(t0, t1, (t1 - t0) / dt + 10)
  u = linspace(t0, t1, (t1 - t0) / dt + 10)
  while r.successful() and r.t < t1:
    r.integrate(r.t + dt)
    e[r.t / dt] = abs(r.y[0][0]) ** 2
    g[r.t / dt] = abs(r.y[1][1]) ** 2
    u[r.t / dt] = abs(r.y[2][2]) ** 2
  return e, g, u


# von Neumann equation's
def right_part(t, rho):
  hamiltonian = (h / 2) * array(
    [[delta, omega_s, omega_p / 2.0 * sin(t * w_p)],
    [omega_s, 0.0, 0.0],
    [omega_p / 2.0 * sin(t * w_p), 0.0, 0.0]],
    dtype=complex128)
  return (dot(hamiltonian, rho) - dot(rho, hamiltonian)) / (1j * h)


def create_integrator():
  r = ode(right_part).set_integrator('zvode', method='bdf', with_jacobian=False)
  psi_init = array([[1.0, 0.0, 0.0],
                   [0.0, 0.0, 0.0],
                   [0.0, 0.0, 0.0]], dtype=complex128)
  t0 = 0
  r.set_initial_value(psi_init, t0)
  return r, t0


def main():
  r, t0 = create_integrator()
  t1 = 10 ** -6
  dt = 10 ** -11
  e, g, u = integrate(r, t0, t1, dt)

main()
4

2 に答える 2

4

このような複雑な行列方程式を処理できるcall のラッパーを作成しましたscipy.integrate.odeint。PYTHON で行列結合微分方程式を解くときに固有値をプロットする方法をodeintw参照してください。行列微分方程式を含む別の問題について。

これは、コードの使用方法を示す簡略化されたバージョンです。(簡単にするために、例からほとんどの定数を削除しました)。

import numpy as np
from odeintw import odeintw


def right_part(rho, t, w_p):
    hamiltonian = (1. / 2) * np.array(
        [[0.1, 0.01, 1.0 / 2.0 * np.sin(t * w_p)],
        [0.01, 0.0, 0.0],
        [1.0 / 2.0 * np.sin(t * w_p), 0.0, 0.0]],
        dtype=np.complex128)
    return (np.dot(hamiltonian, rho) - np.dot(rho, hamiltonian)) / (1j)


psi_init = np.array([[1.0, 0.0, 0.0],
                     [0.0, 0.0, 0.0],
                     [0.0, 0.0, 0.0]], dtype=np.complex128)


t = np.linspace(0, 10, 101)
sol = odeintw(right_part, psi_init, t, args=(0.25,))

sol(101, 3, 3)ソリューションを保持するshape を持つ複雑な numpy 配列になりますrho(t)。最初のインデックスは時間インデックスで、他の 2 つのインデックスは 3x3 行列です。

于 2014-11-04T23:46:02.033 に答える
2

QuTiPには、マスター方程式や Lindblad 減衰項などを使用して、これを行うための優れた積分器がいくつかあります。QuTiP 自体は scipy.odeint の薄いラッパーにすぎませんが、特に優れたドキュメントがあるため、多くのメカニズムがより優れたものになります。

于 2015-05-07T12:44:04.753 に答える