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

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!

CodePudding user response：

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.

CodePudding user response：

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.