Solving simultaneous multivariate polynomial equations with python

Question

edit: the reference I got my equations from contained a couple of errors. I've fixed it here. Solutions might actually make sense now!

When a two layer fluid flows over topography, there exist a number of different solutions depending on the relative size of the flow speed and the wave speed in the fluid.

These are termed 'supercritical', 'subcritical' and 'critical' (the first two I refer to here as 'extra-critical').

The following equations define the bounding lines between critical and extra-critical behaviour in (h, U0) parameter space:

I want to eliminate d_1c (i.e. I don't care what it is) and find solutions to these equations in (h, U_0).

Simplifying factors:

• I only need answers for given d_0
• I do not need exact solutions, just an outline of the solution curves, so this can be solved either analytically or numerically.
• I only want to plot over the region (h, U0) = (0,0) to (0.5, 1).

I'd like to solve this using modules available in the Enthought distribuion (numpy, scipy, sympy), but really don't know where to start. It's the elimination of the variable d1c that really confuses me.

Here are the equations in python:

def eq1(h, U0, d1c, d0=0.1):
    f = (U0) ** 2 * ((d0 ** 2 / d1c ** 3) + (1 - d0) ** 2 / (1 - d1c - d0) ** 3) - 1
    return f

def eq2(h, U0, d1c, d0=0.1):
    f = 0.5 * (U0) ** 2 * ((d0 ** 2 / d1c ** 2) - (1 - d0) ** 2 / (1 - d1c - d0) ** 2) + d1c + (h - d_0)
    return f

I'm expecting a solution that has a number of solution branches (not always physical, but don't worry about that) and looks roughly like this:

How do I go about implementing this?

Answer 1

Semi-formally, the problem you are trying to solve is the following: given d0, solve the logical formula "there exists d1c such that eq1(h, U0, d1c, d0) = eq2(h, U0, d1c, d0) = 0" for h and U0.

There exists an algorithm to reduce the formula to a polynomial equation "P(h, U0) = 0", it's called "quantifier elimination" and it usually relies on another algorithm, "cylindrical algebraic decomposition". Unfortunately, this isn't implemented in sympy (yet).

However, since U0 can easily be eliminated, there are things you can do with sympy to find your answer. Start with

h, U0, d1c, d0 = symbols('h, U0, d1c, d0')
f1 = (U0) ** 2 * ((d0 ** 2 / d1c ** 3) + (1 - d0) ** 2 / (1 - d1c - d0 * h) ** 3) - 1
f2 = U0**2 / 2 * ((d0 ** 2 / d1c ** 2) + (1 - d0) ** 2 / (1 - d1c - d0 * h)) + d1c + d0 * (h - 1)

Now, eliminate U0 from f1 and insert the value in f2 (I'm doing it "by hand" rather than with solve() to get a prettier expression):

U2_val = ((f1 + 1)/U0**2)**-1
f3 = f2.subs(U0**2, U2_val)

f3 only depends on h and d1c. Also, since it's a rational fraction, we only care about when its numerator goes to 0, so we get a single polynomial equation in 2 variables:

p3 = fraction(cancel(f3))

Now, for a given d0, you should be able to invert p3.subs(d0, .1) numerically to get h(d1c), plug it back into U0 and make a parametric plot of (h, U0) as a function of d1c.

Answer 2

Let me deal with elimination of d1c first. Imagine you managed to massage the first equation to have the form d1c = f(U, h, d0). Then you would substitute this into the second equation, and have a certain relation between U, h and d0. With d0 fixed, this defines a single equation for two variables, U and h, from which you can, in principle, find U for any given h. Which seems to be what you are calling solutions, based on your last sketch. Bad news is that it's not easy to get d1c from either of your equations. Good news is that you don't need to.

fsolve can take a system of equations, say two equations which depend on two variables and give you the solution. In this case: fix h (d0 is fixed already), and feed to fsolve the system you have, treating as variables U0 and d1c. Record the value of U0, repeat for next value of h, and so on.

Notice that, contrary to the suggestion of @duffymo, I'd recommend using fsolve, or, at least start with it, and look for other solvers only if it runs out of steam.

A possible caveat is that you are expecting to have more than one solution for U0 given h: fsolve needs a starting guess, and there's no easy way of telling it to converge to one of the solution branches. If that turns out to be a problem, look at brentq solver.

A different way would be to observe that you can easily eliminate U0 from your system. This way you'll get a single equation for h and d1c, solve it for d1c for each value of h, and then use either of your original equations to calculate U0 given d1c and h.

An example of using fsolve:

>>> from scipy.optimize import fsolve
>>> def f(x, p):
...   return x**2 -p
... 
>>> fsolve(f, 0.5, args=(2,))
array([ 1.41421356])
>>> 

Here args=(2,) is the syntax for telling fsolve that what you actually want to solve if f(x,2)=0, and 0.5 is the starting guess for the value of x.