Finding complex roots from many nonlinear equations in python

Question

I tested an algorithm that has been published in the literature, which involves solving a set of "m" nonlinear equations in both Matlab and Python. A set of nonlinear equations includes input variables that contain complex numbers, so the resulting solutions must also be complex. At the moment, I was able to get good results in Matlab 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 the complex matrix and vector, respectively. The solution converges very well for a wide range of initial values.

I tried to imitate these results in Python with very little success. Numerical solvers seem to be configured differently, and the levenburg-marquardt algorithm exists under the root of the function. 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 tried to use some other algorithms, such as "broyden2" and "anderson", but they are much inferior to Matlab and give good results after playing with the initial conditions. The fsolve function also cannot handle complex variables.

I was wondering if there is something that I am applying incorrectly, and if anyone has an idea, maybe how to solve complex nonlinear equations in Python correctly.

Thank you very much

Answer 1

When I come across 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 that takes a complex input array x (let's 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])] 

In fsolve you can use the new function real_f : the real and imaginary parts of the function are simultaneously solved, considering the real and imaginary parts of the input argument as independent.