Numerical Solutions for System of Non-Linear Equation in Python

I have 2 simple equations:

k = 2.013*h^0.4917 and h = 3.57*k^0.4917

These equations can be solved analytically where k = 5.77 and h = 8.47. I tried to solve it in Python using fsolve and I have followed the way from: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html#scipy.optimize.fsolve

Below is my code:

from scipy.optimize import fsolve

def equations(x):
    return[x[0] - 2.013*x[1]**0.4917,x[1] - 3.57*x[0]**0.4917]

root =  fsolve(equations, [1, 1])

And the result is

<ipython-input-133-11ce0ecaa7e4>:2: RuntimeWarning: invalid value encountered in double_scalars
  return[x[0] - 2.013*x[1]**0.4917,x[1] - 3.57*x[0]**0.4917]
RuntimeWarning: The iteration is not making good progress, as measured by the 
  improvement from the last ten iterations.
  warnings.warn(msg, RuntimeWarning)


print(root)
array([1., 1.])

I am not sure what I did wrong here such that I did not get the correct result. Would you like to point my mistake for this problem? Thank you.

Answer

1. fsolve does not know that your variables are non-negative. The solver goes into the negative zone (because from (1, 1) gradients tell to go towards the negative zone), gets NaNs there, and gets stuck. You should tell somehow where you are looking for a solution. fsolve does not support bounds directly. least_squares can do this.

2. There are several solutions. (0, 0) is also fit. You probably have to somehow push away from it, because it seems not what you want to get from a solver.

You may try something like that.

from scipy.optimize import least_squares

def equations(x):
    return[x[0] - 2.013*x[1]**0.4917,x[1] - 3.57*x[0]**0.4917]

root = least_squares(equations, [2, 2], bounds=([1, 1], [10, 10]))
print(root)

x: array([5.74279193, 8.43196966])