menturi628's profile - activity

I am attempting to solve for a solution of a system of two non-linear implicit equations using the following code:

x = var('x')
y = var('y')
P = [(-1,-5), (1,-5), (-5,0), (5,5)]

# Defining the function
d = sum([sqrt( (x-p[0])^2 + (y-p[1])^2 ) for p in P])
show(d)

# Differentiate with respect to x and y
eqx = d.diff(x)
eqy = d.diff(y)

# Plot both implicit curves
g1 = implicit_plot( eqx==0, (x,-10,10), (y,-10,10), color="blue" )
g2 = implicit_plot( eqy==0, (x,-10,10), (y,-10,10), color="red" )
show(g1 + g2) # note that you can clearly see an intersection of the two curves

# Solve for the solution
print("Solving...")
sol = solve([eqx==0, eqy==0], x, y) # this gets stuck or takes a long time
show(sol)

Everything runs, up to the point of the solve function, which continues to run for what appears to be indefinitely. The code show(g1 + g2) shows a graph that clearly shows there exists an intersection for both curves. I tried to use to_poly_solve=True without success. I do not mind an approximate solution, however I was unable to find a numeric solver for a system such as this (find_root afaik only works on one variable) that will work.

Does there exist a numeric solver which is capable of solving a system of this form? What other alternatives are there?

Thanks, menturi

I'm not certain if the behaviour is correct in this example, or I am possibly misunderstanding something. Sorry about the lengthy example, but I had difficulties simplifying it into a simpler one.

# Setup units
lbf  = units.force.pound_force;
inch = units.length.inch;
sec  = units.time.second;
psi  = lbf / inch^2;
gal  = 231.0*inch^3;
min  = 60.0*sec;
gpm  = gal/min;
hole = gpm/sqrt(psi);

# Define the given constraints
FLext = 20000.0*lbf;
vext  = 15.0*inch/sec;
FLret = -10000.0*lbf;
vret  = -15.0*inch/sec;

# Assumed values
Ps = 3000.0*psi;
pv = 1.0;
pc = 1.4;

# System of 3 equations, 3 unknowns
ABE = var('A_BE');
ARE = var('A_RE');
kv = var('k_v');
eq1 = vext^2 == (Ps*ABE - FLext) * kv^2 / ABE^3 / (1 + pv^2 * pc^(-3));
eq2 = vret^2 == (FLret + Ps*ARE) * kv^2 / ARE^3 / (1 + pv^2 * pc^3);
eq3 = pc == ABE / ARE;
sol = solve([eq1, eq2, eq3], [ABE, ARE, kv], solution_dict=True);
sol = sol[1];
show(sol[ABE])
show(sol[ARE])
show(sol[kv])

# Showing the units of various variables
show(sol[kv].convert())
show(hole.convert())
show( (sol[kv]/hole).convert() ) # this should be dimensionless

# Showing the potential bug with unit conversion
show(sol[kv].convert(hole)) # ValueError: Incompatible units

Running the code, the section labelled # Showing the units of various variables, you can see the various units of the relevant variables. The code for show( (sol[kv]/hole).convert() ) should return a dimensionless value, however it still shows units in the result. Looking at the non-numeric portions of it, everything appears to cancel. Running show(sol[kv].convert(hole)) results in ValueError: Incompatible units. I expected a conversion into gpm/sqrt(psi) (i.e. holes).

What concerns me is my calculator (TI-92+) has conversion difficulties with this unit as well, being unable to convert. Is this the expected/normal/ideal behaviour for sage? If so, for what reason(s)? If not, what should I do?

Thank you,

menturi