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
On Thursday, January 23, 2020 at 1:17:36 PM UTC-8, Brent W. Baccala wrote:
>
> 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:
>
> [image: ulu.png]
> 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
>
>
--
You received this message because you are subscribed to the Google Groups
"sage-devel" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/sage-devel/e8dd9584-988a-4400-88e5-605bb5fbb114%40googlegroups.com.
