Solving equation with two square roots



Hi Poul,

El Jueves, 18 de Agosto de 2005 13:40, Poul Riis escribió:
> To be more specific my question is:
>
> Is there a 'stronger' command than 'solve', or can I set some parameter so
> that maxima actually can solve an equation like
> solve(sqrt((x-a*e)^2+y^2)+sqrt((x+a*e)^2+y^2)=2*a,y);

This is a quick solution and I am sure it can be improved in many directions. 
First, define two functions, sqrtxp and isolve, which could be written in a 
separate file, and then call isolve. 

If the equation is not too complicated, isolve only solves for square roots, 
not for cubics or so.

This is a sample session:


-------------  BEGIN MAXIMA -----------------

(%i1) sqrtxp(expr,x):=block([z,w:false],
   if not atom(expr) and op(expr) = "-" then expr: -expr,
   if atom(expr) then return(false),
   if op(expr) = "*"
      then z: args(expr)
      else z: [expr],
   for t in z while not w do
      w: is(not atom(t) and op(t) = 'sqrt and not freeof(x,t)),
   w  )$
(%i2) isolve(eq,x):=block([original,lm,rm,change:true,sol,correct:[]],
   original: eq,
   /* if not an equation, construct it */
   if not op(eq) = "=" then eq: eq=0,

   while change do(
      /* if there are no changes in this loop,        */
      /* then the equation is ready for calling solve */
      change: false,

      /* expand the equation */
      eq: factor(eq),
      eq: expand(eq*denom(lhs(eq))*denom(rhs(eq))),

      /* take left an right terms separately */
      lm: lhs(eq),
      rm: rhs(eq),

      /* scan the left side: if a term has not the          */
      /* unknown under sqrt, it is thrown to the right side */
      if atom(lm)
         then eq: eq - lm
         else (if op(lm) # "+" then lm: [lm],
               for t in args(lm) do
                   if not sqrtxp(t,x)
                      then eq: eq - t
                      else change: true),

      /* now scan the right side: if a term has the unknown */
      /* under sqrt, it is thrown to the left side          */
      if not atom(rm)
         then (if op(rm) # "+" then rm: [rm],
               for t in args(rm) do
                   if sqrtxp(t,x)
                      then (eq: eq - t,
                            change: true)),

      /* Now let's make ^2 */
      eq: expand(eq^2)),

   /* hopefully, the eq. is ready to be solved */
   sol:solve(eq,x),

   /* now, solutions must be checked against input equation */
   for r in sol do
      if lhs(r)=x and freeof(x,rhs(r)) and is(expand(at(original,r)))
         then correct:append(correct,[r]),

   /* The output is a list of two lists: 1st list is the  */
   /* result obtained by solve, and the 2nd list contains */
   /* those solutions which solve the original equation.  */
   /* If you want a more simple output, substitute next   */
   /* command by 'correct' and you'll get legal solutions */
   [sol,correct])$
(%i3) isolve(sqrt(x+1)+sqrt(x+6)=5,x);
(%o3)                         [[x = 3], [x = 3]]
(%i4) isolve(x^2-6*x+9=4*sqrt(x^2-6*x+6),x);
(%o4) [[x = 5, x = 1, x = 3 - 2 sqrt(3), x = 2 sqrt(3) + 3],
                          [x = 5, x = 1, x = 3 - 2 sqrt(3), x = 2 sqrt(3) + 
3]]
(%i5) isolve(2*x-x^2+sqrt(6*x^2-12*x+7)=0,x);
(%o5) [[x = 1 - 2 sqrt(2), x = 2 sqrt(2) + 1, x = 1],
                                        [x = 1 - 2 sqrt(2), x = 2 sqrt(2) + 
1]]
(%i6) isolve(sqrt(r^2-x^2)+sqrt(2*r*x-x^2)=m,x);
                        4      2  2    4     3    2
              m sqrt(3 r  + 2 m  r  - m ) - r  - m  r
(%o6) [[x = - ---------------------------------------,
                               2      2
                            2 r  + 2 m
                                            4      2  2    4     3    2
                                  m sqrt(3 r  + 2 m  r  - m ) + r  + m  r
                              x = ---------------------------------------], 
[]]
                                                   2      2
                                                2 r  + 2 m
(%i7) isolve(sqrt((x-a*e)^2+y^2)+sqrt((x+a*e)^2+y^2)=2*a,y);
                    2  2    2    2  2    2
(%o7) [[y = - sqrt(e  x  - x  - a  e  + a ),
                                                  2  2    2    2  2    2
                                        y = sqrt(e  x  - x  - a  e  + a )], 
[]]

-------------  END MAXIMA -----------------

See the comments at the end of isolve for the interpretation of the result. 
The second list contains the solutions for the 'initial' equation. According 
to that, it seems your equation has no solutions.I have no time to check 
this; it's too late here and tomorrow I fly to Barcelona.

Hope this helps,
Mario.

P.D.: By the way, what's the solution given by your TI89? I'm curious.

-- 
--------------------------
Mario Rodriguez Riotorto
www.biomates.net