Solving simultaneous multivariate polynomial equations with python

math, polynomial-math, python, scipy, sympy

Solution

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.

Problem

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?

Original source