Improving `solve`

[Brent W. Baccala]

I recently had to solve a system of equations, and I was a bit disappointed about how difficult it was, and I wondering how to improve solve to make it easier to use.

Quickly, here's the problem.  It's based on the following geometry:

We're given an ellipse with a thin strip protruding from it.  We want to calculate the parameters of a circle that is tangent to both the plate and the ellipse.  D, d, r1, and r2 are given.  We need to compute x, y, and R.

I set it up easily enough:

sage: var('x,y,R,r1,r2,D,d')
sage: eq1 = (x-D)^2+(y-d/2-R)^2==R^2
sage: eq2 = (x/r1)^2 + (y/r2)^2 ==1
sage: eq3 = (x-D)/(y-d/2-R) == (x/y)*(r2/r1)^2

Equation of the circle, equation of the ellipse, and the tangency condition.

Now what?  None of the following attempts work:

sage: solve([eq1, eq2, eq3, r1==1.1, r2==.7, d==.15, D==2], R)
sage: solve([eq1, eq2, eq3, r1==1.1, r2==.7, d==.15, D==2], [R,x,y])
sage: solve([eq1, eq2, eq3], [R,x,y])
sage: solve([eq1, eq2, eq3], [R,x,y], algorithm='sympy')
sage: solve([eq1, eq2, eq3], [R], algorithm='sympy')
sage: solve([eq1, eq2, eq3, r1==1.1, r2==.7, d==.15, D==2], [R,x,y], algorithm='sympy')
sage: solve([eq1, eq2, eq3, r1==QQ(1.1), r2==QQ(.7), d==QQ(.15), D==2], [R,x,y], algorithm='sympy')
sage: solve([eq1, eq2, eq3, r1==QQ(1.1), r2==QQ(.7), d==QQ(.15), D==2], [R,x,y])

In fact, I couldn't dope out any set of arguments that would get solve to do what I want.

Sage can do the calculation, though, like this:

sage: B.<x,y,R> = QQbar[]
sage: vals = {r1:QQ(1.1), r2:QQ(.7), d:QQ(.15), D:2}
sage: eqs = [B((eq.lhs() - eq.rhs()).numerator().subs(vals)) for eq in [eq1, eq2, eq3]]
sage: I = ideal(eqs)
sage: [r for r in I.variety() if r[R] > 0]

The basic idea is to clear denominators, substitute in for constants as early as possible to simplify the ring, then compute the (hopefully zero-dimensional) variety and apply any inequalities at that point.

The problem is the relative obscurity of these functions.

I'm thinking... should we add a new "algorithm" to solve?  One that will use our polynomial ring code?  Any idea how to come up with a heuristic to select when to prefer that algorithm over Maxima, which I believe is the default?

Sounds like a good project for a junior developer.  Do we have any?

    agape

    brent

[rjf]

1. I suggest you rescale your system by a factor of 100. Then all the

numbers are integers or rationals instead of floats.

2. Maxima solves the system with the parameters as integers in 0.16 seconds on my machine,

at least if I've done it right.

Unfortunately it seems to need floats to do so. However..

there appear to be 12 (complex) solutions, some with tiny imaginary parts.

There seems to be at least one that is nearly real,

R=197.2860635696821,x=1.239194860231853*10^-6*%i+83.16095934479846,y=9.072282470609553*10^-7*%i+45.81930799497882

so re-scaling, it seems R=1.97, x= 0.83, y=0.45.

How does that seem?

Thinking about your suggestion of changing "solve",  it could, of course

be done, but "solve" is such a general term that maybe it should be

broken down into other more specific commands.  e.g.

something involving polynomial-systems and so on.

RJF