[SciPy-User] SCIPY FSOLVE

[Etrade Griffiths]

Hi

I am trying to solve directly a series of equations describing flow 
in a network using FSOLVE but have not had much success so far.  A 
small example is given below.

I have a feeling that the problem may be ill-conditioned because 
there are two sources of small numbers in the equations: one 
associated with small coefficients and the other associated with 
small differences in pressure between adjacent nodes.  The main 
problem seems to be that FSOLVE ends up with estimates of the 
variables that result in raising a negative number to a power, so 
that some of the residuals become -1.#IND in value.  I tried trapping 
this type of error but FSOLVE does not appear to be able to move 
towards the solution.  I would be grateful for any suggestions on how 
to "encourage" FSOLVE.

Thanks in advance

Alun Griffiths

# =========================
#
# Simple model of gathering network
# Uses SCIPY FSOLVE to solve for network pressures directly
#

import math
import scipy.optimize

# Set up the system of network equations

def NetworkEquations(x):

     # Transfer current guess to variables

     p2 = x[0]
     q1 = x[1]
     q2 = x[2]
     q3 = x[3]

     # Define constants

     c = [300.0, 2.00E-5, 1.50E-5]
     n = [0.50, 0.75, 0.70]

     p1 = 200.0
     p3 = 2000.0

     # Define the residuals

     residuals = []

     # ... Kirchoff

     residuals.append( q1 - q2 - q3 )

     # ... pressure drop along pipeline

     curr_err = q1 - (  c[0] * ( p2 ** 2 - p1 ** 2 ) ) ** n[0]
     residuals.append(curr_err)

     # ... pressure drop over well 1

     curr_err = q2 - c[1] * ( p3 ** 2 - p2 ** 2 ) ** n[1]
     residuals.append(curr_err)

     # ... pressure drop over well 2

     curr_err = q3 - c[2] * ( p3 ** 2 - p2 ** 2 ) ** n[2]
     residuals.append(curr_err)

     # All equations evaluated so return residuals

     return residuals


# Solve equations using Broyden's method

x0 = [250.0, 2.0, 1.0, 1.0]
x1 = scipy.optimize.fsolve(NetworkEquations,x0)
print x1