Solution to a system of non-linear equations in R^2

Question

I am trying to find a solution to the following system where f and g are R^2 -> R^2 functions:

f(x1,x2) = (y1,y2)
g(y1,y2) = (x1,x2)

I tried solving it using scipy.optimize.fsolve as follows:

def eqm(vars):
    x1,x2,y1,y2 = vars
    eq1 = f([x1, x2])[0] - y1
    eq2 = f([x1, x2])[1] - y2
    eq3 = g([y1, y2])[0] - x1
    eq4 = g([y1, y2])[1] - x2
    return [eq1, eq2, eq3, eq4]

fsolve(eqm, x0 = [1,0.5,1,0.5])

Although it is returning an output, it does not seem to be a correct one as it does not seem to satisfy the two conditions, and seems to vary a lot with the x0 specified. Also getting a warning: 'The iteration is not making good progress, as measured by the improvement from the last ten iterations.' I do know for a fact that a unique solution exists, which I have obtained algebraically.

Not sure what is going on and if there is a simpler way of solving it, especially using just two equations instead of splitting up into 4. Something like:

def equations(vars):
    X,Y = vars
    eq1 = f(X)-Y
    eq2 = g(Y)-X
    return [eq1, eq2]

fsolve(equations, x0 =[[1,0.5],[1,0.5]])

Suggestions on other modules e.g. sympy are also welcome!

Edit: Following is the specific problem I am trying to solve. The equations are derived from an optimization problem. Basically I am trying to find the intersection of reaction_m and reaction_f such that reaction_m([e_f,tau_f])[1:] = [e_m, tau_m] and reaction_f([e_m,tau_m])[1:] = [e_f, tau_f]

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from scipy.optimize import fsolve

alpha = 0.5
alpha_f = 0.3
alpha_m = 0.3
gamma = 1
kappa_c = 1
w_m = 1
w_f = 0.5
p_e = 1
n=2

def reaction_m(y):
    e_f = y[0]
    tau_f = y[1]
    
    def objective_m(x):
        c_m = x[0]
        e_m = x[1]
        tau_m = x[2]
        return -(np.log(c_m)+ \
            kappa_c*np.log(gamma*tau_m**alpha_m*tau_f**alpha_f*(e_m+e_f)**(1-alpha_f-alpha_m)))
        
    def constraint_m(x):
        c_m = x[0]
        e_m = x[1]
        tau_m = x[2]
        return (1-n*tau_m)*w_m - c_m - n*p_e*e_m 

    x0=[1,1,0.5]

    b = [0,1/n]
    nb = [0, np.inf]
    bnds = [nb,nb,b]
    con = {'type':'eq', 'fun':constraint_m}
    
    sol = minimize(objective_m, x0, method='SLSQP', bounds=bnds, constraints=con)
    return sol.x

def reaction_f(y):
    e_m = y[0]
    tau_m = y[1]
    
    def objective_f(x):
        c_f = x[0]
        e_f = x[1]
        tau_f = x[2]
        return -(np.log(c_f)+ \
            kappa_c*np.log(gamma*tau_m**alpha_m*tau_f**alpha_f*(e_m+e_f)**(1-alpha_f-alpha_m)))
        
    def constraint_f(x):
        c_f = x[0]
        e_f = x[1]
        tau_f = x[2]
        return (1-n*tau_f)*w_f - c_f - n*p_e*e_f 

    x0=[1,1,0.5]

    b = [0,1/n]
    nb = [0, np.inf]
    bnds = [nb,nb,b]
    con = {'type':'eq', 'fun':constraint_f}
    
    sol = minimize(objective_f, x0, method='SLSQP', bounds=bnds, constraints=con)
    return sol.x

def cournot(vars):
    e_m, tau_m, e_f, tau_f = vars
    eq1 = reaction_m([e_f, tau_f])[1] - e_m
    eq2 = reaction_m([e_f, tau_f])[2] - tau_m
    eq3 = reaction_f([e_m, tau_m])[1] - e_f
    eq4 = reaction_f([e_m, tau_m])[2] - tau_f
    return [eq1, eq2, eq3, eq4]

fsolve(cournot, x0 = [1,0.7,0.5,0.2])

Answer 1

First, I recommend working with numpy arrays since manipulating these is simpler than lists.

I've slighlty rewritten your code:

import scipy.optimize as opt

def f(x):
  return x
def g(x):
  return x

def func(vars):
  input = np.array(vars)
  eq1 = f(input[:2]) - input[2:]
  eq2 = g(input[2:]) - input[:2]
  return np.concatenate([eq1, eq2])

root = opt.fsolve(func, [1, 1, 0., 1.2])

print(root)
print(func(root)) # should be close to zeros

What you have should work correctly, so I believe there is something wrong with the equations you're using. If you provide those, I can try to see what may be wrong.

Answer 2

This seems to be more of a problem of numerical mathematics than Python coding. Your functions may have "ugly" behavior around the solution, may be strongly non-linear or contain singularities. We cannot help further without seeing the functions. One thing you might try is to instead solve a system

g(f(x)) - x = 0

and simplify g(f(x)) as much as possible analytically. Then calculate y = f(x) after solving the equation.