how to solve 3 nonlinear equations in python

Question

I have the following system of 3 nonlinear equations that I need to solve:

-xyt + HF = 0

-2xzt + 4yzt - xyt + 4z^2t - M1F = 0

-2xt + 2yt + 4zt - 1 = 0

where x, HF, and M1F are known parameters. Therefore, y,z, and t are the parameters to be calculated.

Attemp to solve the problem:

def equations(p):
    y,z,t = p
    f1 = -x*y*t + HF
    f2 = -2*x*z*t + 4*y*z*t - x*y*t + 4*t*z**2 - M1F
    f3 = -2*x*t + 2*y*t + 4*z*t - 1
    return (f1,f2,f3)

y,z,t = fsolve(equations)

print equations((y,z,t))

But the thing is that if I want to use scipy.optimize.fsolve then I should input an initial guess. In my case, I do not have any initial conditions.

Is there another way to solve 3 nonlinear equations with 3 unknowns in python?

Edit:

It turned out that I have a condition! The condition is that HF > M1F, HF > 0, and M1F > 0.

Answer 1

@Christian, I don't think the equation system can be linearize easily, unlike the post you suggested.

Powell's Hybrid method (optimize.fsolve()) is quite sensitive to initial conditions, so it is very useful if you can come up with a good initial parameter guess. In the following example, we firstly minimize the sum-of-squares of all three equations using Nelder-Mead method (optimize.fmin(), for small problem like OP, this is probably already enough). The resulting parameter vector is then used as the initial guess for optimize.fsolve() to get the final result.

>>> from numpy import *
>>> from scipy import stats
>>> from scipy import optimize
>>> HF, M1F, x=1000.,900.,10.
>>> def f(p):
    return abs(sum(array(equations(p))**2)-0)
>>> optimize.fmin(f, (1.,1.,1.))
Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 131
         Function evaluations: 239
array([ -8.95023217,   9.45274653, -11.1728963 ])
>>> optimize.fsolve(equations, (-8.95023217,   9.45274653, -11.1728963))
array([ -8.95022376,   9.45273632, -11.17290503])
>>> pr=optimize.fsolve(equations, (-8.95023217,   9.45274653, -11.1728963))
>>> equations(pr)
(-7.9580786405131221e-13, -1.2732925824820995e-10, -5.6843418860808015e-14)

The result is pretty good.