Solving system of nonlinear complex equations in Python

Question

I'm trying to solve a problem with 8 unknowns and 8 complex equations. I've tried to use fsolve but I get the error message:

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

From what I've now read fsolve doesn't support complex equations and hence my questions, how would I solve systems of complex non-linear equations in Python?

PS: I've seen the suggestion to split my problem up into imaginary and real part and use fsolve on those separately but that is too cumbersome.

This is the relevant snippet of my code:

A=1

def equations(p):
    B,C,D,F,G,H,I,J = p
    return (
        A+B-C-D,
        1j*k0*(A-B)                                   -1j*k1*(C-D),

        B*exp(1j*k1*a1) + D*exp(-1j*k1*a1)            - F*exp(1j*k0*a1) - G*exp(-1j*k0*a1),
        1j*k1* ( C*exp(1j*k1*a1) - D*exp(-1j*k1*a1) ) - 1j*k0*( F*exp(1j*k0*a1) - G*exp(-1j*k0*a1) ),

        F*exp(1j*k0*a1L) + G*exp(-1j*k0*a1L)          - H*exp(-k2*a1L) - I*exp(k2*a1L),
        1j*k0*( F*exp(1j*k0*a1L) - G*exp(-1j*k0*a1L) )- k2*( -H*exp(-k2*a1L) + I*exp(k2*a1L) ),

        H*exp(-k2*a12L) + I*exp(k2*a12L)              - J*exp(1j*k0*a12L),
        k2*( -H*exp(-k2*a12L) + I*exp(k2*a12L) )      - 1j*k0*J*exp(1j*k0*a12L)
    )

B, C, D, F, G, H, I, J = fsolve(equations, (-1,-1,-1,-1,-1,-1,-1,-1))

`scipy.optimize.fsolve` for finding roots of functions ( https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.optimize.fsolve.html ) — Mattias Johnson, Oct 09 '19 at 10:58

Carllos Limma · Answer 1 · 2022-11-18T14:13:01.620

Algorithm in which the code below was written algorithm for Newton’s Method for Systems


# Author  : Carlos Eduardo da Silva Lima
# Theme   : Newton’s Method for Systems (real or complex)
# Language: Python
# IDE     : Google Colab
# Data    : 18/11/2022

######################################################################
# This Part contains the imports of packages needed for this project #
######################################################################
import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import inv, norm, multi_dot
from scipy.optimize import fsolve

######################
# Enter problem data #
######################
x0  = np.array([1.0+0.0j, 1.0+0.0j, 1.0+0.0j]) # Initial guess for the possible root of the set of nonlinear equations entered in F(x)
TOL = 1e-49                                    # Stipulated minimum tolerance
N   = 10                                       # Number of required maximum iterations

############################
# Newton-Raphson algorithm #
############################
def newtonRapshonSistem(F,J,x0,TOL,N):
  x = x0 # First kick
  k = 1
  while(k<=N):

    # Ccalculation of the product between the inverse of the Jacobian matrix (J(x)^(-1)) and the vector F(x) (T a transpose)
    y = -((inv(J(x)))@(F(x).T))
    x += (y.T)

    # absolute value norm
    erro_abs = np.linalg.norm(np.abs(y))
    if erro_abs<TOL:
      break
    k += 1

  # Exit #
  print(f"Number of iterations: {k}")
  print(f"Absolute error: {erro_abs}\n")
  print("\nSolução\n")
  for l in range(0,np.size(x),1):
    print(f"x[{l}] = {x[l]:.4}\n")
    # print(f'x[{l}] = {np.real(x[l]):.4} + {np.imag(x[l]):.4}j')
  return x
  

#################
# Function F(x) #
#################
def F(x):
  # definition of variables (Arrays)
  x1,x2,x3 = x

  # definition of the set of nonlinear equations
  f1 = x1+x2-10000
  f2 = x1*np.exp(-1j*x3*5) + x2*np.exp(1j*x3*5) - 12000
  f3 = x1*np.exp(-1j*x3*10) + x2*np.exp(1j*x3*10) - 8000

  return np.array([f1, f2, f3], dtype=np.complex128)

############
# Jacobian #
############
def J(x):
  # definition of variables (Arrays)
  x1,x2,x3 = x

  # Jacobean matrix elements
  df1_dx1 = 1
  df1_dx2 = 1
  df1_dx3 = 0

  df2_dx1 = np.exp(-1j*x3*5)
  df2_dx2 = np.exp(1j*x3*5)
  df2_dx3 = x1*(-1j*5)*np.exp(-1j*x3*5)+x2*(1j*5)*np.exp(1j*x3*5)

  df3_dx1 = np.exp(-1j*x3*10)
  df3_dx2 = np.exp(1j*x3*10)
  df3_dx3 = x1*(-1j*10)*np.exp(-1j*x3*10) + x2*(1j*10)*np.exp(1j*x3*10)

  matriz_jacobiana = np.array([
      [df1_dx1, df1_dx2, df1_dx3], 
      [df2_dx1, df2_dx2, df2_dx3], 
      [df3_dx1, df3_dx2, df3_dx3]], dtype=np.complex128)

  return matriz_jacobiana

# Calculate the roots
s = newtonRapshonSistem(F,J,x0,TOL,N)

# Application of the result obtained in x in F.
F(s)

Finally! If you don't agree, or if you find any errors, please let me know. In the most I hope to help you and the community the community. Up until :-).

Solving system of nonlinear complex equations in Python

1 Answers1