help links subscribe changes refresh edit

	FriCAS Library Examples
ExampleIntegration2 ←	↑	→ ExampleSolve2

ExampleSolve

Edit detail for ExampleSolve revision 1 of 1

	1
Editor: test1 Time: 2020/02/10 18:11:11 GMT+0
Note:

changed:
-
This is motivated by discussion in Sage mailing list.  We are given sum of ellipsis and
a strip symmetric with respect to one of the axis of the ellipsis.  We want to
find circle with center on the line x = D which is tangent both to ellipsis and the strip.
More precisely, we want to find common point (x, r) and radius R of the circle.
We get equations (for convenience written as polynomials):
\begin{axiom}
)set output algebra on
)set output tex off
eq1 := (x/r1)^2 + (y/r2)^2 - 1
eq2 := (x-D)^2+(y-d/2-R)^2 - R^2
eq3 := (x-D)/(y-d/2-R) - (x/y)*(r2/r1)^2
\end{axiom}
The first is equation of ellipsis, the second of the circle,
third is tangency condition.

This system can be solved, but it would take long time, the answer is large and would be hard to read.  So first
we replace r, r2, D, d by sensible concrete values:
\begin{axiom}
eq1e := eval(eq1, [r1 = 11/10, r2 = 7/10, d = 15/100, D = 2])
eq2e := eval(eq2, [r1 = 11/10, r2 = 7/10, d = 15/100, D = 2])
eq3e := eval(eq3, [r1 = 11/10, r2 = 7/10, d = 15/100, D = 2])
\end{axiom}
Now we can solve
\begin{axiom}
rese := solve([eq1e, eq2e, eq3e], [x, y, R])
\end{axiom}
Not small, but we see that x and y are expressed in terms of R and
for R we get equation of degree 6.  Let us see numeric solutions
for R:
\begin{axiom}
solve(rese(1)(3), 1.0e-15)
\end{axiom}
So there are two real solutions, positive one solve our problem,
negative represents circle such that has interior (a disc) has intersection
with interior of the strip.  This looks unnatural, but is correct
solution of the equations we gave.  To exclude it we would have
to add an inequality, say y > 0.

This is motivated by discussion in Sage mailing list. We are given sum of ellipsis and a strip symmetric with respect to one of the axis of the ellipsis. We want to find circle with center on the line x = D which is tangent both to ellipsis and the strip. More precisely, we want to find common point (x, r) and radius R of the circle. We get equations (for convenience written as polynomials):

fricas

(1) -> )set output algebra on

fricas

)set output tex off

eq1 := (x/r1)^2 + (y/r2)^2 - 1

          2 2     2 2     2  2
        r1 y  + r2 x  - r1 r2
   (1)  ----------------------
                  2  2
                r1 r2

Type: Fraction(Polynomial(Integer))

fricas

eq2 := (x-D)^2+(y-d/2-R)^2 - R^2

         2                   2           1  2          2
   (2)  y  + (- d - 2 R)y + x  - 2 D x + - d  + R d + D
                                         4

Type: Polynomial(Fraction(Integer))

fricas

eq3 := (x-D)/(y-d/2-R) - (x/y)*(r2/r1)^2

              2     2         2      1         2
        ((- r2  + r1 )x - D r1 )y + (- d + R)r2 x
                                     2
   (3)  -----------------------------------------
                   2 2      1         2
                 r1 y  + (- - d - R)r1 y
                            2

Type: Fraction(Polynomial(Fraction(Integer)))

The first is equation of ellipsis, the second of the circle, third is tangency condition.

This system can be solved, but it would take long time, the answer is large and would be hard to read. So first we replace r, r2, D, d by sensible concrete values:

fricas

eq1e := eval(eq1, [r1 = 11/10, r2 = 7/10, d = 15/100, D = 2])

        100  2   100  2
   (4)  --- y  + --- x  - 1
         49      121

Type: Fraction(Polynomial(Fraction(Integer)))

fricas

eq2e := eval(eq2, [r1 = 11/10, r2 = 7/10, d = 15/100, D = 2])

         2             3      2          3     6409
   (5)  y  + (- 2 R - --)y + x  - 4 x + -- R + ----
                      20                20     1600

Type: Polynomial(Fraction(Integer))

fricas

eq3e := eval(eq3, [r1 = 11/10, r2 = 7/10, d = 15/100, D = 2])

          72             49      147
        (--- x - 2)y + (--- R + ----)x
         121            121     4840
   (6)  ------------------------------
                2           3
               y  + (- R - --)y
                           40

Type: Fraction(Polynomial(Fraction(Integer)))

Now we can solve

fricas

rese := solve([eq1e, eq2e, eq3e], [x, y, R])

   (7)
   [
     [
         x
       = 
             21647170639872527456000   5    563398654925279868893800   4
           -------------------------- R  - -------------------------- R
           11314902455332490895259929      33944707365997472685779787
         + 
            227978377910623361179980  3    543170461396404207065460341  2
           ------------------------- R  - ---------------------------- R
           3771634151777496965086643      1086230635711919125944953184
         + 
            103497888021194504792764569     18719246207620127101109311295
           ---------------------------- R + -----------------------------
           2413845857137598057655451520      7724306742840313784497444864
       ,

         y
       = 
              8578201025986535925690800000   5
           - ------------------------------ R
             481302005742478165211671599873
         + 
               687645572606252172049144000   4
           - ------------------------------ R
             160434001914159388403890533291
         + 
             569404955344370090785647085975  3
           - ------------------------------ R
             962604011484956330423343199746
         + 
           231817205837663178957237643885  2
           ------------------------------ R
           855648010208850071487416177552
         + 
         15047964589394026276064984240695      56533764389483354647866157177923
         -------------------------------- R + ---------------------------------
          6845184081670800571899329420416     273807363266832022875973176816640
       ,

            6   11937  5   16384102559  4   77692853121  3   1498249076221  2
           R  + ----- R  + ----------- R  - ----------- R  - ------------- R
                60760       476358400        4763584000       15243468800
         + 
             66093598324107     2139938711851329
           - -------------- R - ----------------
              3810867200000      48779100160000
       = 
         0
       ]
     ]

Type: List(List(Equation(Fraction(Polynomial(Fraction(Integer))))))

Not small, but we see that x and y are expressed in terms of R and for R we get equation of degree 6. Let us see numeric solutions for R:

fricas

solve(rese(1)(3), 1.0e-15)

   (8)  [R = - 1.4998115943_024088459, R = 1.9728606253_044955743]

Type: List(Equation(Polynomial(Float)))

So there are two real solutions, positive one solve our problem, negative represents circle such that has interior (a disc) has intersection with interior of the strip. This looks unnatural, but is correct solution of the equations we gave. To exclude it we would have to add an inequality, say y > 0.