# Solving a 2nd order complex-valued matrix differential equation in Python

I am trying to solve the following complex-valued matrix differential equation backwards (i.e. not starting at $$r=0$$, but rather at $$r > 0$$):

$$F'' = 2ikF' + VF$$.

Here $$F=F(r)$$ and $$V=V(r)$$ are 2x2 matrices. $$F$$ can have complex entries, $$V$$ cannot.

The boundary conditions are

$$\lim_{r \to \infty} F(r) = \mathbb{1}_2$$

and

$$\lim_{r \to \infty} F'(r) = 0$$.

In practice setting $$F =\mathbb{1}_2$$ at e.g. $$r=10$$ and integrating towards $$r=0$$ should be sufficient to impose this boundary condition.

The problem is that I'm finding it difficult to find a solver that can handle solving complex-valued matrix differential equations backwards - or maybe I just do not understand the interfaces sufficiently.

It seemed that my best bet was to use odeintw, but I am getting the error TypeError: Cannot cast array data from dtype('complex128') to dtype('float64') according to the rule 'safe' error

when defining the function as

def asys(f, r, v, k):
return 2j*k*f + np.matmul(v, f)


An added complication is that $$V$$ is a function of $$r$$, and it is not clear to me how to pass it to the solver.

If it is at all helpful, the equation is related to the scattering theory of two interacting scalar fields with equal masses. I initially solved this equation years ago with hacked together Fortran code that I can no longer find. It was also a nightmare to understand or change anything (it lacked proper abstractions), so it would be great if I could implement a Python solution that is more transparent in its operation ... Can anyone suggest an idea of how to solve this equation using odeintw, or any other Python library?

• Couldn't we convert it to a system of real-valued ODEs? Commented Dec 30, 2022 at 1:30
• @MaximUmansky I would prefer not to. But if you think that is a good solution feel free to post the solution as an answer. Commented Dec 30, 2022 at 9:19
• V(r) is given as data?
– knl
Commented Dec 30, 2022 at 21:53
• @knl No, it is defined as a 2x2 function of r. Typically a constant symmetric real matrix multiplied by something like e**(-mur) or e**(-mur**2). Commented Dec 31, 2022 at 7:41
• I was using some math PDE jargon. Here "data" refers to something that is given by the user and not solved as a part of the solution routine. We like to clearly distinguish between unknowns and data. So I take V is given and F is unknown which wasn't completely clear from your original question.
– knl
Commented Dec 31, 2022 at 16:04

## 1 Answer

This turned out to be embarrassingly simple. I'm posting the answer in the hope that it will help some future person who makes the same mistake.

Basically, the initial condition needs to be complex if the solution is, so you need to explicitly define the initial condition as complex, even if the imaginary part is 0.

Using the following initial condition solved my problem:

    # initial condition
f0 = np.array([
[1.0 + 0j, 0.0],
[0.0, 1.0 + 0j]
])


This 'feature' is listed in the docstring of the odeintw function: https://github.com/WarrenWeckesser/odeintw/blob/master/odeintw/_odeintw.py, but not on the project's main github page, so it's easy to miss if you don't delve into the solver code.