[Bug] NumPy's asarray wrong casting when evoked from SciPy's integrate.solve_ivp

[homocomputeris]

As RK45 supports complex numbers, it should be able to solve a simple Schrödinger equation, right?

It seems that solve_ivp cannot automatically detect than the propagator's type is complex, which will result in a complex vector after the 1st iteration, and then the solver proceeds as if it were still real.

Reproducing code example:

import numpy as np
import scipy as sp

dim = 8
shp = (dim, dim)
# a complex Hermitian matrix
R = np.random.rand(dim, dim) + 1j*np.random.rand(dim, dim)
H = R + R.conj().T
# a real vector
psi = np.random.rand(dim)
# this makes integration possible:
# psi = np.asarray(psi, dtype=np.complex128)


def Hf(t, y):
    return -1j*H@y


sol = sp.integrate.solve_ivp(Hf, (0, 1), psi).y[:, -1]
exact = sp.linalg.expm(-1j*H)@psi

difference = sol - exact

Error message:

/usr/lib/python3.6/site-packages/numpy/core/numeric.py:492: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)

Scipy/Numpy/Python version information:

1.0.0 1.14.0 sys.version_info(major=3, minor=6, micro=4, releaselevel='final', serial=0)

[homocomputeris]

Trace's interesting part:

  File "/home/nik/Documents/Code/GniPy/scipy_bug.py", line 36, in <module>
    sol = sp.integrate.solve_ivp(Hf, (0, 1), psi).y[:, -1]
  File "/usr/lib/python3.6/site-packages/scipy/integrate/_ivp/ivp.py", line 455, in solve_ivp
    solver = method(fun, t0, y0, tf, vectorized=vectorized, **options)
  File "/usr/lib/python3.6/site-packages/scipy/integrate/_ivp/rk.py", line 99, in __init__
    self.f = self.fun(self.t, self.y)
  File "/usr/lib/python3.6/site-packages/scipy/integrate/_ivp/base.py", line 139, in fun
    return self.fun_single(t, y)
  File "/usr/lib/python3.6/site-packages/scipy/integrate/_ivp/base.py", line 21, in fun_wrapped
    return np.asarray(fun(t, y), dtype=dtype)
  File "/usr/lib/python3.6/site-packages/numpy/core/numeric.py", line 492, in asarray
    return array(a, dtype, copy=False, order=order)
  File "/usr/lib/python3.6/warnings.py", line 99, in _showwarnmsg
    msg.file, msg.line)
  File "/home/nik/Documents/Code/GniPy/krypy_bug.py", line 18, in warn_with_traceback
    traceback.print_stack(file=log)
/usr/lib/python3.6/site-packages/numpy/core/numeric.py:492: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)

[pv]

Whether the equation is considered complex-valued is probably determined by the data type of `psi`, instead of the solver trying to guess it based on return values of the callable. The documentation states " For problems in a complex domain pass `y0` with a complex data type (even if the initial guess is purely real)."

[homocomputeris]

Wouldn't it be better to use numpy.common_type(y0, H(0, y0)) or whatever? What are the problems it may cause?
I think automatic type definition can be expected. We write programs to offload tasks from the brain, not to create more.

[drhagen]

The feature to support complex-valued solvers was proposed in the development repository, but never implemented.

[derdotte]

I have recently come across this problem aswell. This time with the dirac equation. This should definitely get a revisit and perhaps the current status quo should be reevaluated.

EDIT: 7 December:
A friend of mine told me just a few hours ago that he had a problem with wrong casting aswell while solving his system (this time it was a free basis expansion of the schrödinger equation for calculating hyperfinestructure transition probabilities). This really needs a reevaluation. Might need a better error message if the behavior isnt going to change, might need to be more prominant in the documentation. Eitherway this is definitely not how a python library should operate, after all a python library user should mostly not need to care about inputing the correct combination of types such that the expected operation of the function is achieved.

[mdhaber]

Just thought I'd mention a relevant discussion on the mailing list.