[SciPy-User] SCIPY FSOLVE

Warren Weckesser warren.weckesser@enthought....
Mon Oct 12 08:50:12 CDT 2009


It appears you are correct in your concern about the closeness of p1 and 
p2.  You can avoid the problem of raising a negative number to a 
fractional power if you reformulate your residual.

If you change this:

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

to this:

     curr_err = q1**(1/n[0]) - (  c[0] * ( p2 ** 2 - p1 ** 2 ) )

the script gives the answer:

[ 200.00004794    2.39840655    1.77542113    0.62298543]

Note that p2 is very close to p1, so it is not surprising that during 
the solver's iterations, p2**2 - p1**2 occasionally becomes negative.

Warren


Etrade Griffiths wrote:
> 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
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>   



More information about the SciPy-User mailing list