Solve a non-linear system of equations

Question

I am trying to solve this non-linear system of equations:

from scipy.optimize import fsolve
import numpy as np

def equations(vars):
    x, y, z = vars
    eq1 = x - 0.095*np.log(z) - 1.2022
    eq2 = y - 100000/x
    eq3 = z - y
    return [eq1, eq2, eq3]

ini = (1, 1, 1)
x, y, z =  fsolve(equations, ini)

print(x, y, z)

The system gives me a solution but is not the solution to the problem. The solver gives me this solution: x=2.22015, y=14373.01967, z=14372.9181 but the real solution is x=2.220157, y=45041.83986, z=45041.83986.

It seems that the problem is the initialization of the values. If I put these values for initialization:

ini = (2, 40000, 40000)
x, y, z =  fsolve(equations, in)

The system gives me the real solution: x=2.220157, y=45041.83986, z=45041.83986

What can I do to obtain the good solution without knowing it in advance?

Answer 1

Try this, it loops thru 3 ranges for ini, call solve and if status is 1 we return because status 1 is a success or pass status. We set full_output parameter to true in fsolve() to get status info.

Code

import time
from scipy.optimize import fsolve
import numpy as np


def equations(vars):
    x, y, z = vars
    eq1 = x - 0.095*np.log(z) - 1.2022
    eq2 = y - 100000/x
    eq3 = z - y
    return [eq1, eq2, eq3]


def sol():
    ret = None
    for i in range(1, 1000):
        print(f'processing ... {i}')

        for j in range(1, 1000):
            for k in range(1, 1000):
                ini = (i, j, k)
                res =  fsolve(equations, ini, full_output=True)
                if res[2] == 1:  # if status is 1 then it is solved.
                    return res
                ret = res

    return ret

# Test
t1 = time.perf_counter()

res = sol()
print(f'status: {res[2]}, sol: {res[0]}')
print(f'elapse: {time.perf_counter() - t1:0.1f}s')

Output

processing ... 1
status: 1, sol: [2.22015798e+00 4.50418399e+04 4.50418399e+04]
elapse: 2.9s