Numerical simulations are never perfect. Of course, the fundamental reason for this is the limit set by floating-point numbers. However, more often than not, it is an error made by the user in coding their problem that leads to the primary source of error in a calculation. Here, we will look at one such error by considering the dynamics of a driven system composed of two-coupled oscillators. Here we drive oscillator B with a classical drive, and assume the initial state is both oscillators in the ground state. The QuTiP code, and resulting plot is given below:
"""
Driven system of two-coupled harmonic oscillators
"""
from qutip import *
from pylab import *
wa = 1.0 * 2 * pi # frequency of system a
wb = 1.0 * 2 * pi # frequency of system a
wab = 0.2 * 2 * pi # coupling frequency
ga = 0.05 * 2 * pi # dissipation rate of system a
gb = 0.05 * 2 * pi # dissipation rate of system b
Na = 7 # number of states in system a
Nb = 7 # number of states in system b
E= 1 * 2 * pi
a = tensor(destroy(Na), qeye(Nb))
b = tensor(qeye(Na), destroy(Nb))
na = a.dag() * a
nb = b.dag() * b
H = wa*na+wb*nb+wab*(a.dag()*b+a*b.dag())+E*(b+b.dag())
# start with both oscillators in ground state
psi0 = tensor(basis(Na), basis(Nb))
#construct collapse operators
c_op_list = []
c_op_list.append(sqrt(ga) * a)
c_op_list.append(sqrt(gb) * b)
tlist = linspace(0, 5, 101)
#run simulation
data1 = mesolve(H,psi0,tlist,c_op_list,[na,nb])
#plot results
plot(tlist,data1.expect[0],'b',tlist,data1.expect[1],'r',lw=2)
xlabel('Time',fontsize=14)
ylabel('Excitations',fontsize=14)
legend(('Oscillator A', 'Oscillator B'))
show()
Fig1: Coupled oscillator evolution.
Now, looking at the plot, we see that the results are what we would expect intuitively from our initial setup. Since we drive oscillator B, we would expect the number of excitations in B to grow large before the coupling with oscillator A starts to drive energy into this mode. At later times, we see that the energy goes back and forth.
If this graph was in a publication, chances are it would pass the referees just fine. Especially if we are like everyone else and do not attach the source code for this plot! However, the graph above is in fact wrong. Our model is correct, so qualitatively it gives reasonable answers. However, quantitatively it is off at some times by over 100%! So where did we screw up? In our initial setup, we chose to have N=7 states for both of our resonators. If we were naive, then looking at the plot, we might assume that this is a big enough Hilbert space since we have at most ~3 excitations in any given oscillator. Fortunately for us, we paid attention in class, and we remember that an harmonic oscillator driven by a classical drive tends to produce coherent states, and a coherent state has a Poisson distribution over a wide range of number states. Therefore, in choosing only 7 states, we are in fact chopping off some of the system dynamics. Since this truncation will be most pronounced when the system has the largest number of excitations. lets look at what a more accurate model with N=10 states gives us at time t=0.50 for oscillator B:
 |
Fig.2: Distribution for oscillator B at t=0.50
overlaid with coherent state distribution.
Blue states are truncated when N=7.
|
From Fig. 2, we can immediately see the error. We are effectively ignoring the contributions from the higher modes (blue) of the generated (nearly) coherent state when N=7. It should not be surprising that the results at later times are not correct. Including more states (N=10) gives the following results.
 |
| Fig 3.: Comparison of dynamics with N=7 and N=10 states. |
In this case, the average error is 22.4%, and as alluded to earlier, is over 100% at some time-steps. We can look at the error between the data for N states vs. N-1 as we increase N:
 |
| Fig. 4: Error between data for N and N-1 oscillator states. |
We can see that the error begins to level off starting around N=10, where the increase in the number of oscillator states has little effect on the numerics.
The final message: Check to see if the number of states you are using is in fact large enough by visually checking the effect of adding more states, or better yet, plotting the error as a function of system size. And as always, include your code in your publications.
