The instances of `+` (to be used for arithmetic) called within the procedures of the LinearAlgebra package proper are almost all instances involving the binding to the global :-`+`. Those bindings occur at the time of creation and saving of the procedures to the Maple library.

You would likely need to redefine global :-`+` at least, in order to get a custom procedure for `+` to be used within LinearAlgebra:-MatrixMatrixMultiply. And even that might not be enough to succeed, if the work is passed off to some kernel builtin like mvMultiply.

An overloaded rebinding of `+` to your new procedure is not by itself going to suffice for what you seem to be describing.

How about overloading your own `.`, to use your `+` for matrix-matrix multiplication?

(I don't follow your description of your `+` for your one brief, incomplete example. Also, what about linearity?)

How about,

   frontend(collect, [four, [V^P, V^(P+1), V^(P+2)]],
            [{`*`, `+`, list}, {}]);

or,

   frontend(collect, [four, [seq(V^(P+i),i=1..5)]],
            [{`*`, `+`, list}, {}]);

Do you need code to figure out the upper bound of that `seq`?

Following others ideas below. The solve command needs quite a bit of hand-holding here.

restart:
expr:=sin(Pi*(x+1)/(4*x^2-4*x+2))=cos(Pi*(x-2)/(4*x^2-4*x+2)):
new:=simplify(trigsubs(convert((lhs-rhs)(expr),sin)),size)[];

                /        2     \    /   / 2        \\
                |    Pi x      |    |Pi \x  - x - 1/|
          -2 cos|--------------| sin|---------------|
                |   2          |    |   2           |
                \4 x  - 4 x + 2/    \4 x  - 4 x + 2 /

solsA:={solve({new,x<2},x,allsolutions,explicit)}:
solsB:={solve({new,x<2},x)}:
solsC:={solve({new,x>=2},x)}:
maybes:=solsA union solsB union solsC:

sols,others:=selectremove(type,maybes,set(identical(x)=realcons)):

sols;

    /          /    1   1   (1/2)\    /    1   1   (1/2)\   
   { {x = 1}, { x = - - -- 5      }, { x = - + -- 5      }, 
    \          \    2   10       /    \    2   10       /   

      /      1  (1/2)   1\    /    1  (1/2)   1\ \ 
     { x = - - 5      + - }, { x = - 5      + - } }
      \      2          2/    \    2          2/ / 

# How best to show that x::real in others[1] and others[2] only when _Z1=0 (x=1)?
others[1];

            /                                 (1/2)\ 
            |                /      2        \     | 
           <     2 _Z1 + 1 + \-4 _Z1  - 2 _Z1/      >
            |x = ----------------------------------| 
            \                4 _Z1 + 1             / 

seq(about(t), t in indets(others[1],name) minus {x,constants});

Originally _Z1, renamed _Z1~:
  is assumed to be: integer

others[2];

          /                                    (1/2)\ 
          |                   /      2        \     | 
         <       -2 _Z1 - 1 + \-4 _Z1  - 2 _Z1/      >
          |x = - -----------------------------------| 
          \                   4 _Z1 + 1             / 

seq(about(t), t in indets(others[2],name) minus {x,constants});

Originally _Z1, renamed _Z1~:
  is assumed to be: AndProp(integer,OrProp(RealRange(-1/2,Open(-1/4)),RealRange(Open(-1/10),0)))

Have you already tried resizing the inlined plot window manually? Ie, click-select it with the mouse pointer, and then drag its highlighted outer edges to resize it.

Alternatively... how are you exporting the plot? By right-click and context-menu on the plot? You may be able to get by by using the `plotsetup` command to export the plot programmatically. This provides options for the height and width.

restart:
with(plots):

p1:=plot(x^2,x=0..1,scaling=constrained,view=[-4..5,-1..3]):

plotsetup(ps, plotoutput=cat(kernelopts(homedir),
                             "/My Documents/bypd.ps"),
          plotoptions="noborder,height=200,width=500");

p1; # this will not be displayed in the worksheet, getting sent to the file instead

plotsetup(default);

p1; # this will be displayed in the worksheet

Of course the following will not work on all possible examples. But it works here.

solve({x^2+y^2+z^2 = 3, x+y+z = 3}, {x,y,z}):

q:=[allvalues(%)]:

{seq(solve(evalc(`union`((Im~,Re~)(Q)))),Q=q)}[];

                            {x = 1, y = 1, z = 1}

For multivariable problems fsolve uses Newton's method, after checking a couple of things such as whether the system is linear.

interface(verboseproc=3):
eval(`fsolve/gensys`);
eval(`fsolve/genroot`);
eval(`fsolve/sysnewton`);

It could be strengthened for underdetermined systems. If fsolve were to get support for multivariable constrained root-finding then it would likely have to rely upon more sophisticated optimization techniques, as the DirectSearch package does.

Here is one way.

remdup:=proc(L::list)
  local i,T;
  T:=table([]);
  for i from 1 to nops(L) do
    if not member(L[i],T) then
      T[i]:=L[i];
    end if;
  end do;
  [entries(T,nolist)]; 
end proc:

The above procedure's do-loop can also be done using `seq` or `map`. For example, the do-loop could be replaced by either of

seq(`if`(not member(e,T),assign(T[e],e),NULL), e in L);

map(e->`if`(not member(e,T),assign(T[e],e),NULL),L);

But those don't seem to perform that much better.

If anyone has something that uses much less resources it'd be interesting to see. I tested with

M,N:=5000,5000:
L:=convert(LinearAlgebra:-RandomVector(N,generator=1..M),list):

I suppose that some techniques that ListTools:-MakeUnique uses might be offlimits too. Its following method is very efficient. It is similar to the `map` replacement for the do-loop shown above. But it does faster table lookups and NULL assignments, doing them on the remember table of procedure `once` and not making overt calls to `assign` and `member`.

proc(L::list)
  local once;
  once:=proc(x) once(x):=NULL; x; end proc:
  map(once,L);
end proc:

method=_NCrule may have enough accuracy (even after 5 iterations), and be faster for the embedded numerical integral.

Download CEV_idea.mw

proc(n::posint)
  local i;
  add( (-1)^(i-1) * 1/i^2 , i=1..n );
end proc:

There are many ways to do this. Another simple way, without using the add() command, is

proc(n::posint)
  local i, result;
  result := 0;
  for i from 1 to n do
    result := result + (-1)^(i-1) * 1/i^2;
  end do;
  return result;
end proc:

If "pickup" means the time duration until the helicopter will be able to lift off, and if the "fastest pickup" means the shortest such duration, then the pickup is determined by the longest (max) time taken amongst all four groups. That's as opposed to the sum of the times taken for all four groups, as per your solution.

So, assuming as you did that they all travel the same constant speed,

> max(evalf([distance(a, x), distance(b, x), distance(c, x), distance(d, x)]));

                          19.92133289

> max(evalf([distance(a, o), distance(b, o), distance(c, o), distance(d, o)]));

                          18.38477631

The helicopter would then be able to take off sooner if they all met at the centroid than if they met at the intersection of the two cross-diagonals.

But, in this interpretation, how about meeting at approximately x=13.97282, y=7.054256 for least maximal-amongst-the-four distance travelled of 15.65264 km?

Basically, I think that this sentence is not right, "Granted one group will have to travel farther than the rest but as a whole the groups need to regroup at the intersection for minimal total travelled distance." They should regroup as a whole in the shortest amount of time, and that is determined by the max of the four distances travelled.

You could try constructing and testing the combinations one at a time.

> with(combinat):                                

> A:=[a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u]:
> M,N := nops(A), 10:                            

> B := [$1..N]:                                  
> testthis:=A[[op(B)]];                          

                  testthis := [a, b, c, d, e, f, g, h, i, j]

> for ii from 1 to 3 do ### to numbcomb(M, N)-1 do   
>    B:=nextcomb(B, M);                          
>    testthis:=A[[op(B)]];                       
> end do;                                        

                     B := {1, 2, 3, 4, 5, 6, 7, 8, 9, 11}

                  testthis := [a, b, c, d, e, f, g, h, i, k]

                     B := {1, 2, 3, 4, 5, 6, 7, 8, 9, 12}

                  testthis := [a, b, c, d, e, f, g, h, i, l]

                     B := {1, 2, 3, 4, 5, 6, 7, 8, 9, 13}

                  testthis := [a, b, c, d, e, f, g, h, i, m]

You might want to use a full colon on the `end do`, to suppress printed output, if you do the full set above.

From the current version of the Programming Manual (section 7.4, subsection "The evalf Command"):

      All floating-point computations are performed in finite precision, with intermediate
      results generally being rounded to Digits precision. As such, it is possible for
      round-off errors to accumulate in long computations. Maple ensures that the result
      of any single floating-point arithmetic operation (+, -, *, /, or sqrt) is fully
      accurate.  Further, many of the basic functions in Maple, such as the trigonometric
      functions and their inverses, the exponential and logarithmic functions, and some of
      the other standard special functions for mathematics, are accurate to within .6 units
      in the last place (ulps)...

Yours is not a single floating-point computation (a call to exp, sqrt, `+`, etc) but is rather a composition of operations. The above paragraph describes that round-off error can accumulate. This has been Maple's general model of floating-point computation for at least the past fifteen years.

Try reducing the ranges allowed for the free variables (once W is fixed), in the call to fsolve. Otherwise when fsolve fails to find a root then Q will appear to be discontinuous. And Optimization does much better for smooth functions.

By reducing the ranges there is a better chance that a randomly generated starting point of the multivariable Newton's method solver might succeed. This can become more key as the number of variables gets larger. Of course, it only works if there is a solution inside the hyperbox.

fsolve(eval(sys,w=W),{a1=0..1,a2=0..0.1,a3=0..1,a4=0..1,a5=0..1,a6=0..1,t=0..1});

With these ranges I got an even plot of Q and an objective value much smaller (near zero, for w near 1.25).

Download maple_qq.mw

You can form this using the Matrix constructor itself along with a procedure that acts according to a conditional.

> restart:
                                
> B:=Matrix([[2,0,0],[0,8,0],[0,0,4]]):   

> Matrix(op(1,B),(i,j)->`if`(i=j,1.5^B[i,j],0));

                        [2.25         0           0   ]
                        [                             ]
                        [ 0      25.62890625      0   ]
                        [                             ]
                        [ 0           0         5.0625]

For your diagonal Matrix B that result would be the same as from

Matrix(op(1,B),(i,j)->`if`(i=j,1.5^B[i,j],B[i,j]));

since the fallback value of B[i,j] will be zero for entries with i not equal to j.

Rather than wrestle with cos and sin terms that decline to simplify easily (shifts, etc, arising in results returned from `solve`), you could generate the results yourself,

restart:

z := -sqrt(6)-sqrt(18)*I;

                        (1/2)        (1/2)
                      -6      - 3 I 2    
 
r,th := abs(z), simplify(arctan(evalc(Im(z)),evalc(Re(z))));

                           (1/2)    2   
                        2 6     , - - Pi
                                    3   

ans := [seq(surd(r,3)*(cos+I*sin)((th+2*k*Pi)/3),k=0..3-1)];

          [ (1/3)  (1/6) /   /2   \        /2   \\  
          [2      6      |cos|- Pi| - I sin|- Pi||, 
          [              \   \9   /        \9   //  

             (1/3)  (1/6) /   /4   \        /4   \\  
            2      6      |cos|- Pi| + I sin|- Pi||, 
                          \   \9   /        \9   //  

             (1/3)  (1/6) /    /1   \        /1   \\]
            2      6      |-cos|- Pi| - I sin|- Pi||]
                          \    \9   /        \9   //]

{seq(simplify(p^3), p in ans)};

                 /  (1/2)  (1/2)        (1/2)\ 
                { -2      3      - 3 I 2      }
                 \                           / 

simplify(map(convert,ans,polar));

  [     / (1/2)  (1/6)    2   \       / (1/2)  (1/6)  4   \  
  [polar|2      3     , - - Pi|, polar|2      3     , - Pi|, 
  [     \                 9   /       \               9   /  

         / (1/2)  (1/6)    8   \]
    polar|2      3     , - - Pi|]
         \                 9   /]

This comes up quite often, that the nth roots of a (exact rational) complex number get computed by `solve` as a mess. It even can even take a few moments to compute them. Couldn't that special form of request be recognized by `solve` and handled immediately and cleanly by just returning the simple sequence? Wouldn't it help many students?