[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