Issues with complex numbers when solving ode in python

Question

I have the following code:

import numpy as np
from matplotlib import pyplot as plt
from scipy.integrate import odeint

exp = np.exp
pi = np.pi

def func(y, t, W, omega, delta, omega_Rabi):
    c_up, c_down = y
    dydt = [-1j*np.exp(-1j*delta*t)*(1j*W+omega_Rabi)*c_down,-1j*np.exp(1j*delta*t)*(-1j*W+omega_Rabi)*c_up]
    return dydt

W = 2*pi*10
omega_Rabi = 2*pi*300;
delta = 2*pi*5000;
omega = 2*pi*100000;
y0 = [0,1]

t = np.linspace(0, 0.01, 100)
sol = odeint(func, y0, t, args=(W,omega,delta,omega_Rabi))

plt.plot(t, abs(sol[:, 0])**2, 'b')
plt.show()

The code runs and gives me a plot, but it is the wrong one. I am getting this warning:

ComplexWarning: Casting complex values to real discards the imaginary part
  output = _odepack.odeint(func, y0, t, args, Dfun, col_deriv, ml, mu,

so it seems to have something to do with the fact that I am using complex numbers, but I don't know what is it or how to fix it. Can someone help me? Thank you!

That is not possible with `odeint`. Of the old solvers you had to use `ode` with method `zvode` to use complex state vectors. Or simply use the newer `solve_ivp` with any method, in the newest version it should work correctly (there was some problem with the step size mechanism). — Lutz Lehmann, Oct 15 '22 at 03:50
Or explicitly make the equation real. Have the state vector with twice the dimension, at the start of the ODE function combine pairs of reals into complex numbers, at the end split those again as pairs of real numbers in a flat array. It may be easier to have the first half of the real array as real parts, the second as imaginary parts. `z=u[:N]+1j*u[N:]` and `du = concatenate([dz.real, dz.imag])` — Lutz Lehmann, Oct 15 '22 at 03:59
@LutzLehmann solve_ivp seems to not give that warning, but the error accumulates so fast. Even if I give a fine spacing for t_eval doesn't change much. — JohnDoe122, Oct 15 '22 at 04:12
`t_eval` does not influence the step size, for that use `atol, rtol` and `max_step`. /// With frequency 5000 you want on the interval to 0.01 at least 200, better 300 or 500 points in a plot to get all parts of the oscillation in full amplitude independent of the phase. — Lutz Lehmann, Oct 15 '22 at 04:24

score 0 · Answer 1 · answered Oct 15 '22 at 03:04

You could just extract the real portion of the complex number by using the .real operation. Following is a revised version of your function.

def func(y, t, W, omega, delta, omega_Rabi):
    c_up, c_down = y
    dydt = [(-1j*np.exp(-1j*delta*t)*(1j*W+omega_Rabi)*c_down).real, (-1j*np.exp(1j*delta*t)*(-1j*W+omega_Rabi)*c_up).real]
    return dydt

Note the use of the .real extraction for each element.

Running your program with that revision still produced the same plot without the warning.

Here is a link you might check for further information. Complex Numbers

Give that a try.

But I don't want just the real part (that is the plot I get, too and it is wrong). I need the imaginary part, too. — JohnDoe122, Oct 15 '22 at 03:56

Issues with complex numbers when solving ode in python

1 Answers1