# QuTiP VS RK45: Which one gives the correct results for time-dependent systems?

I am writing a code for a quantum thermal machine which includes both coherent and dissipative time evolutions in its different stages of operation. However, evolving the system with "mesolve" function of QuTiP (which uses either Adams or BDF methods) and Runge-Kutta gives different results. My code has time dependent Hamiltonians and collapse operators. For demonstration purposes I wrote a simple code which deals with the master equation of a time-dependent harmonic oscillator of mass and frequency equal to 1. The jump operator is the position operator times a time-dependent coefficient.

import numpy as np
import scipy as sp
from qutip import *

tend = 5
tvec = np.linspace(0,tend,1000)
ω1 = 1 #initial frequency
ω2 = 2 #final frequency
gamma = 1 #dissipation rate

def omg_fun(t, args): #time-dependent frequency
s = t/tend
ω_t = ω1 + 10*(ω2-ω1)*s**3 - 15*(ω2-ω1)*s**4 + 6*(ω2-ω1)*s**5
return ω_t

def col_fun(t, args): #time dependent collapse operator coefficient
return np.sqrt(gamma * np.exp(-t))

def omg_fun_rk(t): #time-dependent frequency for RK45
s = t/tend
ω_t = ω1 + 10*(ω2-ω1)*s**3 - 15*(ω2-ω1)*s**4 + 6*(ω2-ω1)*s**5
return ω_t

def col_fun_rk(t):  #time dependent collapse operator coefficient for RK45
return np.sqrt(gamma * np.exp(-t))

N=10
m = 1 #mass
ħ = 1
x = np.sqrt(ħ/(2*m*ω1))*(create(N)+destroy(N)) #position operator
p = 1j*np.sqrt((ħ*m*ω1)/2)*(create(N)-destroy(N)) #momentum operator
H = [(p**2)/(2*m), [(1/2)*m*x**2,omg_fun]]
c_ops = [[x,col_fun]]
rho_init = fock_dm(N,1)
rhot_vec = mesolve(H,rho_init,tvec,c_ops)
rhot_qutip = rhot_vec.states[-1]

#Right hand side of the master equation
def rhs(ρ,H,gamma,op):
c_op = np.sqrt(gamma)*op
res = -1j*(H*ρ - ρ*H)+(1/2)*(2*c_op*ρ*c_op.dag()-ρ*c_op.dag()*c_op-c_op*c_op.dag()*ρ)
return res

#Runge-Kutta Solver
def rksolver(ρ,dt,H,gamma,op):
k1 = rhs(ρ,H,gamma,op)
k2 = rhs(ρ+(dt/2)*k1,H,gamma,op)
k3 = rhs(ρ+(dt/2)*k2,H,gamma,op)
k4 = rhs(ρ+dt*k3,H,gamma,op)
rho = ρ+dt*(k1+2*k2+2*k3+k4)/6
return rho

dt = tvec[1]-tvec[0]
for i in range(len(tvec)):
H_rk = (p**2)/(2*m) + (1/2)*m*x**2*omg_fun_rk(tvec[i])
gamma_rk = col_fun_rk(tvec[i])
rhot_rk = rksolver(rho_init,dt,H_rk,gamma_rk,x)
rho_init = rhot_rk


The final density matrix using mesolve is:

$$\begin{equation*}\left(\begin{array}{*{11}c}0.212 & 0.0 & (0.004-0.052j) & 0.0 & (0.006+0.022j) & 0.0 & (-0.005-0.005j) & 0.0 & (0.002+9.285\times10^{-04}j) & 0.0\\0.0 & 0.384 & 0.0 & (-0.070+0.018j) & 0.0 & (0.009-0.004j) & 0.0 & (0.002+0.001j) & 0.0 & (-0.005-0.002j)\\(0.004+0.052j) & 0.0 & 0.195 & 0.0 & (-0.054+0.023j) & 0.0 & (0.011-0.010j) & 0.0 & (-0.001+0.003j) & 0.0\\0.0 & (-0.070-0.018j) & 0.0 & 0.082 & 0.0 & (-0.039+0.011j) & 0.0 & (0.015-0.008j) & 0.0 & (-0.012+0.011j)\\(0.006-0.022j) & 0.0 & (-0.054-0.023j) & 0.0 & 0.049 & 0.0 & (-0.019+0.009j) & 0.0 & (0.007-0.005j) & 0.0\\0.0 & (0.009+0.004j) & 0.0 & (-0.039-0.011j) & 0.0 & 0.039 & 0.0 & (-0.013-0.002j) & 0.0 & (0.008-0.004j)\\(-0.005+0.005j) & 0.0 & (0.011+0.010j) & 0.0 & (-0.019-0.009j) & 0.0 & 0.019 & 0.0 & (-0.006-0.002j) & 0.0\\0.0 & (0.002-0.001j) & 0.0 & (0.015+0.008j) & 0.0 & (-0.013+0.002j) & 0.0 & 0.008 & 0.0 & (-0.004+5.540\times10^{-04}j)\\(0.002-9.285\times10^{-04}j) & 0.0 & (-0.001-0.003j) & 0.0 & (0.007+0.005j) & 0.0 & (-0.006+0.002j) & 0.0 & 0.005 & 0.0\\0.0 & (-0.005+0.002j) & 0.0 & (-0.012-0.011j) & 0.0 & (0.008+0.004j) & 0.0 & (-0.004-5.540\times10^{-04}j) & 0.0 & 0.007\\\end{array}\right)\end{equation*}$$ The final density matrix using RK45 is:

$$\begin{equation*}\left(\begin{array}{*{11}c}0.243 & 0.0 & (-0.010-0.017j) & 0.0 & (0.003+0.006j) & 0.0 & (-0.001-0.001j) & 0.0 & (3.398\times10^{-04}+5.923\times10^{-05}j) & 0.0\\0.0 & 0.271 & 0.0 & (-0.038+0.003j) & 0.0 & (0.004+0.002j) & 0.0 & (0.002-0.002j) & 0.0 & (-0.004+0.003j)\\(-0.010+0.017j) & 0.0 & 0.186 & 0.0 & (-0.040+0.009j) & 0.0 & (0.008-0.003j) & 0.0 & (-0.001-1.722\times10^{-06}j) & 0.0\\0.0 & (-0.038-0.003j) & 0.0 & 0.105 & 0.0 & (-0.036+0.009j) & 0.0 & (0.013-0.008j) & 0.0 & (-0.013+0.010j)\\(0.003-0.006j) & 0.0 & (-0.040-0.009j) & 0.0 & 0.066 & 0.0 & (-0.020+0.006j) & 0.0 & (0.005-0.005j) & 0.0\\0.0 & (0.004-0.002j) & 0.0 & (-0.036-0.009j) & 0.0 & 0.046 & 0.0 & (-0.014-4.931\times10^{-04}j) & 0.0 & (0.002-0.006j)\\(-0.001+0.001j) & 0.0 & (0.008+0.003j) & 0.0 & (-0.020-0.006j) & 0.0 & 0.030 & 0.0 & (-0.010-0.003j) & 0.0\\0.0 & (0.002+0.002j) & 0.0 & (0.013+0.008j) & 0.0 & (-0.014+4.931\times10^{-04}j) & 0.0 & 0.016 & 0.0 & (-0.005+0.001j)\\(3.398\times10^{-04}-5.923\times10^{-05}j) & 0.0 & (-0.001+1.722\times10^{-06}j) & 0.0 & (0.005+0.005j) & 0.0 & (-0.010+0.003j) & 0.0 & 0.014 & 0.0\\0.0 & (-0.004-0.003j) & 0.0 & (-0.013-0.010j) & 0.0 & (0.002+0.006j) & 0.0 & (-0.005-0.001j) & 0.0 & 0.022\\\end{array}\right)\end{equation*}$$

The fidelity of the final density matrices for these two methods is 0.9778358587429837. Basically they are not close enough. So the question is, which one is correct?!

Note: Changing the solver options and subdividing the time into smaller intervals doesn't change anything. For example, here are the results for the same code but with the following time vector and mesolve options:

tvec = np.linspace(0,tend,10000)
options = Options(nsteps=5000)


The result using mesolve is:

$$\begin{equation*}\left(\begin{array}{*{11}c}0.212 & 0.0 & (0.004-0.052j) & 0.0 & (0.006+0.022j) & 0.0 & (-0.005-0.005j) & 0.0 & (0.002+9.285\times10^{-04}j) & 0.0\\0.0 & 0.384 & 0.0 & (-0.070+0.018j) & 0.0 & (0.009-0.004j) & 0.0 & (0.002+0.001j) & 0.0 & (-0.005-0.002j)\\(0.004+0.052j) & 0.0 & 0.195 & 0.0 & (-0.054+0.023j) & 0.0 & (0.011-0.010j) & 0.0 & (-0.001+0.003j) & 0.0\\0.0 & (-0.070-0.018j) & 0.0 & 0.082 & 0.0 & (-0.039+0.011j) & 0.0 & (0.015-0.008j) & 0.0 & (-0.012+0.011j)\\(0.006-0.022j) & 0.0 & (-0.054-0.023j) & 0.0 & 0.049 & 0.0 & (-0.019+0.009j) & 0.0 & (0.007-0.005j) & 0.0\\0.0 & (0.009+0.004j) & 0.0 & (-0.039-0.011j) & 0.0 & 0.039 & 0.0 & (-0.013-0.002j) & 0.0 & (0.008-0.004j)\\(-0.005+0.005j) & 0.0 & (0.011+0.010j) & 0.0 & (-0.019-0.009j) & 0.0 & 0.019 & 0.0 & (-0.006-0.002j) & 0.0\\0.0 & (0.002-0.001j) & 0.0 & (0.015+0.008j) & 0.0 & (-0.013+0.002j) & 0.0 & 0.008 & 0.0 & (-0.004+5.540\times10^{-04}j)\\(0.002-9.285\times10^{-04}j) & 0.0 & (-0.001-0.003j) & 0.0 & (0.007+0.005j) & 0.0 & (-0.006+0.002j) & 0.0 & 0.005 & 0.0\\0.0 & (-0.005+0.002j) & 0.0 & (-0.012-0.011j) & 0.0 & (0.008+0.004j) & 0.0 & (-0.004-5.540\times10^{-04}j) & 0.0 & 0.007\\\end{array}\right)\end{equation*}$$

The result using RK45 is:

$$\begin{equation*}\left(\begin{array}{*{11}c}0.243 & 0.0 & (-0.010-0.017j) & 0.0 & (0.003+0.006j) & 0.0 & (-0.001-0.001j) & 0.0 & (3.368\times10^{-04}+4.977\times10^{-05}j) & 0.0\\0.0 & 0.271 & 0.0 & (-0.038+0.003j) & 0.0 & (0.004+0.002j) & 0.0 & (0.002-0.002j) & 0.0 & (-0.004+0.003j)\\(-0.010+0.017j) & 0.0 & 0.186 & 0.0 & (-0.040+0.009j) & 0.0 & (0.007-0.003j) & 0.0 & (-0.001+4.236\times10^{-06}j) & 0.0\\0.0 & (-0.038-0.003j) & 0.0 & 0.105 & 0.0 & (-0.036+0.009j) & 0.0 & (0.013-0.008j) & 0.0 & (-0.013+0.010j)\\(0.003-0.006j) & 0.0 & (-0.040-0.009j) & 0.0 & 0.066 & 0.0 & (-0.020+0.006j) & 0.0 & (0.005-0.005j) & 0.0\\0.0 & (0.004-0.002j) & 0.0 & (-0.036-0.009j) & 0.0 & 0.046 & 0.0 & (-0.014-4.878\times10^{-04}j) & 0.0 & (0.002-0.006j)\\(-0.001+0.001j) & 0.0 & (0.007+0.003j) & 0.0 & (-0.020-0.006j) & 0.0 & 0.030 & 0.0 & (-0.010-0.003j) & 0.0\\0.0 & (0.002+0.002j) & 0.0 & (0.013+0.008j) & 0.0 & (-0.014+4.878\times10^{-04}j) & 0.0 & 0.016 & 0.0 & (-0.005+0.001j)\\(3.368\times10^{-04}-4.977\times10^{-05}j) & 0.0 & (-0.001-4.236\times10^{-06}j) & 0.0 & (0.005+0.005j) & 0.0 & (-0.010+0.003j) & 0.0 & 0.014 & 0.0\\0.0 & (-0.004-0.003j) & 0.0 & (-0.013-0.010j) & 0.0 & (0.002+0.006j) & 0.0 & (-0.005-0.001j) & 0.0 & 0.022\\\end{array}\right)\end{equation*}$$

Their fidelity is 0.9778975808571535. So there is no convergence issue.

• If you are interested, I also asked the same question on QuTiP's GitHub. The link is here. github.com/qutip/qutip/issues/1733 Dec 6 '21 at 11:22