Good morning everyone,

I am trying to solve a system of non-linear equations using python (under-expanded jet with losses – Molkov for interested physicists)

The system of equation is shown on the picture below:

In order to do that, I am trying to use the library Scipy and the module fsolve as shown below.

def equations(p):

    P2, T2, u2, rho2, P3, T3, u3, rho3 = p;

    return (P2 - TANK_PRES + rho2 * (u2**2) * ((K / 4) + 1),
            Cp_H2 * TANK_TEMP - Cp_H2 * T2 - ((K + 1) * (u2**2) / 2),
            P2 - ((rho2 * R_H2 * T2) / (1 - ABEL_NOBLE_COEF * rho2)),
            P3 - P2 + rho2 * (u2**2) * ((F / 4) - 1) + rho3 * (u3**2) * ((F / 4) + 1),
            Cp_H2 * T2 + ((u2**2) / 2) - Cp_H2 * T3 - ((F / 4) + 1) * ((u3**2) / 2),
            rho3 - (P3 / (ABEL_NOBLE_COEF * P3 + R_H2 * T3)),
            rho2 * u2 - rho3 * u3,
            u3 - np.sqrt((H2_GAMMA * R_H2 * T3)) / (1 - ABEL_NOBLE_COEF * rho3))


P2, T2, u2, rho2, P3, T3, u3, rho3 =  fsolve(equations, (10**7, 1000, 1500, 1, 10**7, 1000, 1500, 1))
print(equations((P2, T2, u2, rho2, P3, T3, u3, rho3)))
print("*****")
print(P2, T2, u2, rho2, P3, T3, u3, rho3)

But, I encounter a problem due to the np.sqrt((H2_GAMMA * R_H2 * T3)) part in the last equation. Python is trying to solve the system trying a negative value for T3 (and therefore a mathematic error due to negative term with the square-root).
Below the error message received

(-27468750.0, -1967406.7266, -31562027.65265894, 55125.0, -13781.25, 0.7579822540902079, 0.0, -6157.533689399985)
*****
10000000.0 10000.0 1500.0 1.0 10000000.0 10000.0 1500.0 1.0
c:usersguillaume.thirietdownloadstest_h2jet.py:89: RuntimeWarning: invalid value encountered in sqrt
  u3 - np.sqrt((H2_GAMMA * R_H2 * T3)) / (1 - ABEL_NOBLE_COEF * rho3))
c:usersguillaume.thirietdownloadstest_h2jet.py:92: RuntimeWarning: The iteration is not making good progress, as measured by the 
  improvement from the last ten iterations.

I don’t know how to solve the problem. The idea would be to be able to set-up boundaries (min / max values) but it seems fsolve doesn’t allow to do that (or I don’t find anything about that).

I found on the Internet that it is maybe possible to use the least-squares function of Scipy but it does not give me an appropriate result (maybe I am also using it badly).

Can someone help me to solve this system please or provide some tips to solve my issues (other functions from others libraries maybe….)?

Thank you by advance for your help