Cannot plot without high machine precision problems on small numbers

[aphobic]

Hi, I am trying to plot the solution of the Lindblad equation, for a Bell state probability density matrix, however I get non-zero entropy at t=0, which has to be zero.

All calculations show that the y-axis value of the plot should be 0 at t=0, but I still get a plot that starts at 0.032.

I think this is caused by numerical precision problems, of very very small numbers, such as 1.11022302e-16 . Why this is nonzero is strange it self, however, how can I round-up this to give zero ?

The code is attached



Thanks

[deanhystad]

I think the problem is here:

# Evolve the systems and compute entropy at each time step
for t in times:
    rho1 = lindblad_step(rho1, dt, gamma, H, C)
    entropy_values_psi_plus.append(entropy(rho1))
    
    rho2 = lindblad_step(rho2, dt, gamma, H, C)
    entropy_values_psi_minus.append(entropy(rho2))

This does not give rho1 or rho2 values for t == 0. The first rho1/2 value is for t == dt.

There's also a problem here:

# Time evolution
time_steps = 1000
dt = 0.01
times = np.linspace(0, time_steps * dt, time_steps)

I think you want times to be [0, 0.01, 0.02, ... 9.99, 10.0], but what you get is [0, 0.01001001, 0.02002002...9.97997998, 9.98998999, 10.0].

It is difficult to use a for loop because of how you compute rho. I might write your loop like this:

dt = 0.01
rho1 = rho_psi_plus
rho2 = rho_psi_minus
times = []
entropy_values_psi_plus = []
entropy_values_psi_minus = []
for t in range(1001):
    times.append(t * dt)
    entropy_values_psi_plus.append(entropy(rho1))
    entropy_values_psi_minus.append(entropy(rho2))
    rho1 = lindblad_step(rho1, dt, gamma, H, C)
    rho2 = lindblad_step(rho2, dt, gamma, H, C)

Or this:

dt = 0.01
rho1 = rho_psi_plus
rho2 = rho_psi_minus
times = np.linspace(0.0, 10.0, 1001)
entropy_values_psi_plus = []
entropy_values_psi_minus = []
for _ in times:
    entropy_values_psi_plus.append(entropy(rho1))
    entropy_values_psi_minus.append(entropy(rho2))
    rho1 = lindblad_step(rho1, dt, gamma, H, C)
    rho2 = lindblad_step(rho2, dt, gamma, H, C)

The important thing is if you want to plot a point at time == 0, you need to compute a point for time == 0 instead of starting at time = dt.