Finding complex roots from set of non-linear equations in python

Question

I have been testing an algorithm that has been published in literature that involves solving a set of 'm' non-linear equations in both Matlab and Python. The set of non-linear equations involves input variables that contain complex numbers, and therefore the resulting solutions should also be complex. As of now, I have been able to get pretty good results in Matlab by using the following lines of code:

lambdas0 = ones(1,m)*1e-5;
options = optimset('Algorithm','levenberg-marquardt',...
'MaxFunEvals',1000000,'MaxIter',10000,'TolX',1e-20,...
'TolFun',1e-20);

Eq = @(lambda)maxentfun(lambda,m,h,g);
[lambdasf]  = fsolve(Eq,lambdas0,options);

where h and g are a complex matrix and vector, respectively. The solution converges very well for a wide range of initial values.

I have been trying to mimic these results in Python with very little success however. The numerical solvers seem to be set up much differently, and the 'levenburg-marquardt' algorithm exists under the function root. In python this algorithm cannot handle complex roots, and when I run the following lines:

lambdas0 = np.ones(m)*1e-5

sol = root(maxentfun, lambdas0, args = (m,h,g), method='lm', tol = 1e-20, options = {'maxiter':10000, 'xtol':1e-20})

lambdasf = sol.x

I get the following error:

minpack.error: Result from function call is not a proper array of floats.

I have tried using some of the other algorithms, such as 'broyden2' and 'anderson', but they are much inferior to Matlab, and only give okay results after playing around with the initial conditions. The function 'fsolve' also cannot handle complex variables either.

I was wondering if there is something I am applying incorrectly, and if anybody has an idea on maybe how to properly solve complex non-linear equations in Python.

Thank you very much

Answer 1

When I encounter this type of problem I try to rewrite my function as an array of real and imaginary parts. For example, if f is your function which takes complex input array x (say x has size 2, for simplicity)

from numpy import *
def f(x):
    # Takes a complex-valued vector of size 2 and outputs a complex-valued vector of size 2
    return [x[0]-3*x[1]+1j+2, x[0]+x[1]]  # <-- for example

def real_f(x1):
    # converts a real-valued vector of size 4 to a complex-valued vector of size 2
    # outputs a real-valued vector of size 4
    x = [x1[0]+1j*x1[1],x1[2]+1j*x1[3]]
    actual_f = f(x)
    return [real(actual_f[0]),imag(actual_f[0]),real(actual_f[1]),imag(actual_f[1])]

The new function, real_f can be used in fsolve: the real and imaginary parts of the function are simultaneously solved for, treating the real and imaginary parts of the input argument as independent.

Answer 2

Here append() and extend() methods can be used to make it automatic and easily extendable to N number of variables

def real_eqns(y1):
y=[]
for i in range(N):
    y.append(y1[2*i+0]+1j*y1[2*i+1])
real_eqns1 = eqns(y)
real_eqns=[]
for i in range(N):
    real_eqns.extend([real_eqns1[i].real,real_eqns1[i].imag])
return real_eqns