How to interpret a solution from Maple's $solve$ command when I use $assume$?

Question

Consider the following function

$$W(k,a,b,q_1,q_2,q_3,q_4)=\frac{k}{a}(q_1q_3+q_2q_4)+\frac{k}{b}(q_1q_2+q_3q_4)+\frac{k}{\sqrt{a^2+b^2}}(q_2q_3+q_1q_4)$$

In Maple,

W := (k, a, b, q__1, q__2, q__3, q__4) -> k*(q__1*q__3 + q__2*q__4)/a + k*(q__1*q__2 + q__3*q__4)/b + k*(q__2*q__3 + q__1*q__4)/sqrt(a^2 + b^2)

For context, if $q_i$ is an electric point charge, $a$ and $b$ are the sides of a rectangle, and $k$ is Coulomb's constant, then $W$ represents the work required to bring the four charges from an infinite distance away from each other to the corners of the rectangle.

Assume that two charges are negative and two are positive, and the charges have same magnitude $c>0$.

As an example, say $q_1=q_2=c>0$, and $q_3=q_4=-c$

$$W(k,a,b,e,e,-e,-e)=-\frac{2ke^2}{a}+\frac{2ke^2}{b}-\frac{2ke^2}{\sqrt{a^2+b^2}}$$

$$=2ke^2 \left ( \frac{1}{b}-\frac{1}{a}-\frac{1}{\sqrt{a^2+b^2}} \right )$$

I'd like to know if this work can be positive.

By simple inspection it seems there are solutions. If $a$ is very large the $\frac{1}{a}$ and $\frac{1}{\sqrt{a^2+b^2}}$ are very small, so only the positive term $\frac{1}{b}$ remains.

In Maple,

assume(0 < k, 0 < a, 0 < b, 0 < e)
W2 := W(k, a, b, e, e, -e, -e)
solve(0 < W2)

Which results in

{b = b, e = e, k = k, a < 0}

This seems to indicate that there is not solution, because I assumed that $a>0$. But given that I told Maple to assume this, how come it gave me this answer?

Answer 1

It's not clear what kind of solution you are expecting. Note that if there are infinitely many solutions (in the same dimension as the number of variables) then the best you should expect from solve is a terser representation of the region. In that situation solve won't provide you with sample/witness points.

Consider your expression W2.

W:=(k,a,b,e) -> -2*k*e^2/a + 2*k*e^2/b
                -2*k*e^2/sqrt(a^2+b^2):
W2 := collect(W(k,a,b,e), [k,e]);
    /  2   2         2       \  2  

W2 := |- - + - - --------------| e  k
        |  a   b            (1/2)|

        |          / 2    2\     |

        \          \a  + b /     /

For the above form of W2, if we know that e>0 and k>0 then one terser representation would involve the term,

W3 := eval(W2, [e=1,k=1]);
         2   2         2       
 W3 := - - + - - --------------
         a   b            (1/2)
                 / 2    2\     
                 \a  + b /

since it follows that would also have to be positive.

Let's look at a plot:

plots:-implicitplot(W3>0,a=0..1,b=0..1,
    axes=boxed,, size=[500,300],
    filledregions,scaling=constrained);

If the upper edge of that region has a "nicer" representation than W3>0 for that W3 expression then that's about as good as we might expect. So, let's see what we can figure out.

I'm going to deduce a new formula with that radical term isolated on one side.

eq := R = W3:
P := indets(eq,radical)[1];
                  1       
       P := --------------
                      (1/2)
              / 2    2\     
              \b  + a /     


eq2 := P = thaw(solve(eval(eq,P=freeze(P)),
                      freeze(P)));
           1            R a b - 2 a + 2 b

eq2 := -------------- = - -----------------
                  (1/2)           2 a b

         / 2    2\

         \b  + a /

Now I will square both sides. That gets rid of the radical, though I must be careful later because solving the ensuing new equation may produce "solutions" which do not satisfy the original (due to a sign difference).

eq3 := map(`^`,eq2,2);
                                  2
        1      (R a b - 2 a + 2 b) 

eq3 := ------- = --------------------
          2    2            2  2

         b  + a          4 a  b

Now I have a new system to solve, including inequalities on the unknown names. I'll remove any R=... from the result, because I already know R in terms of a and b. There are several solutions from this call to solve, and I'll take the first.

foo := remove(type,
              solve({eq3,R>0,a>0,b>0},
                    [a,b,R])[1],
              identical(R)=anything);
foo := [0 < a,
                 /  4       3       2  2         3    4                 \
       b &lt; RootOf\_Z  - 2 _Z  a + _Z  a  - 2 _Z a  + a , index = real[1]/,
       0 &lt; b]

That form of so-called "indexed RootOf" is not common. It stands for the "first" real root of that quartic. It seems that Maple didn't find a nice explicit representation for it, even given the knowledge that parameter a>0. Perhaps we can do better.

zp := op(1,indets(foo,RootOf)[1]);
      4       3       2  2         3    4

zp := _Z  - 2 _Z  a + _Z  a  - 2 _Z a  + a
S := [solve(zp,real,explicit)];
   [ /     /1   1  (1/2)   1   \         \   

S := [{ _Z = |- + - 2      - - %1| a, a = a }, 
       [ \     \2   2          2   /         /
 /     /1   1  (1/2)   1   \         \ ]     
{ _Z = |- + - 2      + - %1| a, a = a }]     
 \     \2   2          2   /         / ]     


                            (1/2)                     
             /        (1/2)\                          
       %1 := \-1 + 2 2     /

The root we're interested in is this one:

eval(_Z,S[1]);
/                                (1/2)\

   |1   1  (1/2)   1 /        (1/2)\     |

   |- + - 2      - - -1 + 2 2     /     | a
   \2   2          2                     /

Great. That's the root of _Z in the above RootOf. How nice, it's linear w.r.t parameter a.

We can substitute back into the solution foo.

ans := subs(indets(foo,RootOf)[1]=eval(_Z,S[1]),
            foo);
     [           /                                (1/2)\         ]
     [           |1   1  (1/2)   1 /        (1/2)\     |         ]

ans := [0 < a, b < |- + - 2      - - -1 + 2 2     /     | a, 0 < b]
         [           \2   2          2                     /         ]

And that's it, the nicer/terser representation of the shaded region in that plot, which is the solution space we wanted.

We might hope that the is facility could check the result. But it does not (Maple 2022.0)

is(W3>0) assuming op(ans); # weakness
        FAIL

It can check this subregion with another linear upper boundary that doesn't have that radical form,

is(eval(W3,a=2*b)>0) assuming a>0, b>0;
        true

Lastly, if you enter the following into www.wolframalpha.com then you can produce the same terser solution:

solve({-2/a + 2/b - 2/sqrt(b^2 + a^2)>0, a>0, b>0})

Answer 2

solve does not interact well with assume. Instead, give the constraints to solve directly (even with lots of options trying to kick it into the right place):

W := (k, a, b, q__1, q__2, q__3, q__4) -> k*(q__1*q__3 + q__2*q__4)/a + k*(q__1*q__2 + q__3*q__4)/b + k*(q__2*q__3 + q__1*q__4)/sqrt(a^2 + b^2):
constraints := 0 < k, 0 < a, 0 < b, 0 < e:
W2 := W(k, a, b, e, e, -e, -e):
solve({0<W2, constraints},{b,e,k,a}, allsolutions=true,conditionalsolutions=true,explicit=false,tryhard=true);

and Maple warns "solutions may have been lost" while returning no solutions.