Go to ?procedure and read about the uses clause of a procedure header. You can't use with inside a procedure, and I'd recommend against writing procedures that rely on a with performed outside the procedure. With the uses clause, you can make abbreviations:

proc()
   uses St= Statistics, LA= LinearAlgebra;
   local A,B, ...;
   ...;
   A:= St:-Rank(...);
   B:= LA:-Rank(...);
   ...
end proc;

Also, my opinion is that A:-B looks more elegant than A[B]. Maple already has enough square brackets; it is nice to balance the variety of symbols.

(A prior version of this answer was messed up due to the well-known bug of the MaplePrimes editor dropping angle brackets.)

The below assumes that the factors are represented as 2x4 Matrices, and it represents the product as a 4x4 Matrix.

Download tensor.mw

Check out ?ToInert and ?dismantle. The former returns a result that is logically identical to what you want; it is just not presented visually as a tree. The latter displays the output in tree form, but includes details that you might find irrelevant. It would be relatively easy to convert the ToInert output to an explicitly presented tree by using the GraphTheory package and its DrawNetwork command.

Example:

ToInert(x+2*y*z^2);

Excellent questions, TomM. There is no bug, but the answers to your questions are very poorly documented, if documented at all. It took me several hours on 17-Jan-2013 to figure out the answers to those same questions.

TomM wrote:

On Maple 16.02, plattools[transform] only works for very simple procedures. It does not work for procedures  with the `if` function even with the simple arrow definition. or with any if-block in a more general proc definition.

You can use procedures of arbitrary complexity, but it requires a little trick at the beginning of the procedure. Let's call the procedure which is passed to plottools:-transform "the point transformer", and let's call the procedure that plottools:-transform returns "the plot transformer".  When the plot transformer is called, the first call that it makes to the point transformer will be with symbolic (i.e. non-numeric) arguments. It does this to check the number of dimensions returned (because plottools:-transform can be used to transform 2D to 3D and vice versa). On this first call, the point transformer must return a list (or other object) with the correct number of elements.

TomM wrote:

simp:= plot(sin(x), x= 0..Pi);
f:=(xx,yy)->`if`(yy > 0, [xx,4*yy], [xx,yy]);
trns:=plottools:-transform(f);
trns(simp);

So, for the above, your point transformer should be written:

f:=(xx,yy)-> `if`(not (xx::numeric and yy::numeric), [xx,yy], `if`(yy > 0, [xx,4*yy], [xx,yy]));

For an example of a much more complex point transformer, see my MaplePrimes post BubblePlot plotting view: Modification so that the axes don't change.

TomM wrote:

From examining op(trns), the error seems to occur in the procedure: `transform/object`. However, I cannot find a way to print out this procedure, even with verboseproc=2, with print or op. Where is this procedure, and how can I print it out?

It is a procedure local to module plottools. So, to print it out, do

kernelopts(opaquemodules= false);  #Get access to module locals
interface(verboseproc= 2);
eval(plottools:-`transform/object`);

You can also showstat that procedure, or trace it; which is how I figured out the answer to your first problem.

TomM wrote:

On the other hand, with f:=(xx,yy)-> [xx, 3*max(0,yy)+yy], the transform works and produces the same plot as the non-working sequence above. 

It works because the above f will return a list of two elements even if it is passed symbolic arguments. But if it was inside an `if`, then it would just return the unevaluated `if`, which is not a list, and doesn't have two elements.

P:= n-> plot([x^2,2*x-1], x= 1-.1^n..1+.1^n):
plots:-display(<P(0)|P(1)|P(2)>);

You can take some first baby steps towards debugging thus:

kernelopts(opaquemodules= false);
interface(verboseproc= 3);
trace(RegularChains:-TRDisolate_real_zeros);
showstat(RegularChains:-TRDisolate_real_zeros);
printlevel:= 2;

And then run your command. But it gets deep rather quick: RegularChains has 1381 local procedures---the most I've ever seen---all of whose names begin TRD.

For lists G and K, it can be done with elementwise operators thus:

Ga := Basis({(a*~G)[], ((1-a)*~K)[]}, 'tord', deglex(a,r,u,v,w));

The problem is that you have more unknowns than equations. You should look at the output instead of ending every command with a colon. Using N=3, M=3, I get 12 "U" variables and only 6 equations in Eqs1. I think you might need to incorporate some BCs into the equations, but I can't help you much with that. You also might have a problem with index numbers.

@Preben Alsholm There are 32 unknowns and 24 equations. The number of solutions might depend on which 8 variables solve decides to use in the basis, and that may change on different runs.

I tried to solve for arbitrary matrix T as the Asker's updated post asks. I used eliminate. I let it run for 2 hours and use 7 Gigs memory, and then I killed it.

Try raising your cutoff value of 2.5e-2 to 0.3 or 0.4. Then you can really see the points from which it is interpolating the surface. You can see a collection of points on each side of the surface such that the surface lies neatly between them. It might also be instructive to vary the symbolsize of your plotted points, with points closer to the surface plotted larger. This might give more of a clue as to how the algorithm works.

The trick is to make each frame include a static (i.e. not insequence) display of all previous frames. This can be done with any animation. Here's the original code (which I indented) and my modifications. I added a splash of color. I also changed the sum to add for efficiency, which necessitated making g a procedure.

(***** begin original (comment out) *****
g:= sum(arctan(1/sqrt(n)), n = 1 .. i):
p:= seq(
    plots[polygonplot](
        [[0, 0],
        [sqrt(k+1)*cos(subs(i = k, g)), sqrt(k+1)*sin(subs(i = k, g))],
        [sqrt(k)*cos(subs(i = k-1, g)), sqrt(k)*sin(subs(i = k-1, g))]
        ], color = white
    ), k = 1 .. 70
):
plots[display]([p], insequence = true);
***** end original (comment out) *****)

N:= 70:
g:= k-> add(arctan(1/sqrt(n)), n= 1..k):
p:= [seq](
    plots:-polygonplot(
        [[0, 0], sqrt(k+1)*~[cos,sin](g(k)), sqrt(k)*~[cos,sin](g(k-1))    
        ], color= COLOUR(HUE, k/N)
    ), k= 1..N
):
plots:-display([seq](plots:-display(p[1..K]), K= 1..N), insequence = true);

The trick to using evalhf here is that the integration needs to be done at the time h is created, not at the time that h is called. And there's extra complexity that n = 2 needs to be handled as a separate case to avoid division by zero. So, you need to define h this way:

h := subs(J= int(f(x)*sin(n*x), x= -Pi..Pi)/Pi, n-> `if`(n=2, 1, J));

Your f and g can remain the same. If you define h as above, then the call evalhf(g(2)) will work. This is correct usage, not effective usage. To use evalhf effectively, you need to pack as much computation as possible into each invocation of evalhf.

I can't answer the "why?" of your question, but there is a very easy workaround: Make the order 2 or more.

Unfortunately, freeze doesn't work for integers and rationals. But subs does work. So the following kludge does it:

ex:= 16^j/16;
exf:= subs(_x= 16, simplify(subs([1/16= 1/_x, 16= _x], ex)));

It sure looks ugly though.

You have a line

H := dsolve*(ivp, {e(t), w(t)}, numeric);

The asterisk in that line should be removed.

Also, it is not an error, but there is no good reason to use `union` in prefix form. You might as well write

ivp:= sys union ic;

which is closer to the usual mathematical notation.

There are many ways to plot it. Here are two:

plots:-odeplot(H, [[t,e(t)],[t,w(t)]], t= 0..100, color= [blue,red]);

plots:-odeplot(H, [e(t),w(t)], t= 0..100);