[SciPy-User] SCIPY FSOLVE
josef.pktd@gmai...
josef.pktd@gmai...
Mon Oct 12 08:12:16 CDT 2009
On Mon, Oct 12, 2009 at 7:55 AM, Etrade Griffiths
<etrade.griffiths@dsl.pipex.com> 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
>
I never used optimize.fsolve, but in your case, I would try
to reparameterize p2 so it has to be between between
p1 and p3,
p1<=p2<=p3
maybe as a fraction in interval [0,1] and then transform it
to the interval [p1,p3]
There might be other ways to impose the constraint, but
my first attempt is usually reparameterization.
Josef
More information about the SciPy-User
mailing list