@J4James Your original plot3d call produces a GRID data substructure within the PLOT3D data structure itself.

On the other hand, when the plot3d call is modified to use a variable bound for the inner range then the result is a PLOT3D structure containing a MESH.

Ie, considering just one of the two surfaces from your example,

restart:
with(plots):
lambda1:=(S+sqrt(S^2+4*alpha))/(2):
P0:=plot3d(lambda1, alpha=-5..5, S=-10..10):
P1:=plot3d(lambda1, alpha=max(-1/4*S^2,-5)..5, S=-10..10):
op(P0);
op(P1);

It is the GRID vs MESH aspect which is making the difference to the nature of the resulting surface's own grid. Sorry, but I don't know of any source which discusses the rendering of these structures in great detail.

My call to `isolate` was just a way to get a parametrization of the joining spacecure between the two surface pieces. It was fortune, in a sense, that formulaically it worked out simply for this example. I wanted the (formulaic) variable end-point of the range just because I happened to recall how the patching/gridding of surfaces could get done, according to whether it is GRID or MESH.

Another example is below, (where I use some thought to figure out which result from `solve` might suffice, though perhaps the mimimum under assumptions of realness might simplify to what was wanted here). I flip the x- and y-axes, just to get the same exact view as the original.

F12 := -(1/2*(S+sqrt(S^2+4*alpha)))*alpha:
F22 := (1/2*(-S+sqrt(S^2+4*alpha)))*alpha:
solve(F12=F22,alpha);
plot3d({F12,F22}, S=-20..20, alpha=-1..0);
plottools:-transform((x,y,z)->[y,x,z])(
   plot3d({F12,F22}, alpha=max(-1,-S^2/4)..0, S=-20..20));

I'd be surprised if there wasn't an even easier way to force a suitable MESH in the plot3d result.

@Axel Vogt Yes, I was using 64bit Maple 17.02 on Windows 7.

I find it interesting that evalb returns true for the "equality" test of two procedures for which addressof returns two different results.

@John Fredsted 

For HermitianTranspose of Matrix A, there is A^%H

@madssdam The `with` command is not supposed to be utilized within a procedure. Consider utilizing use or `uses` instead (the latter of which is mentioned in the proc help-page) for adjusting the environment within the scope of that procedure.

Are you hoping to change the global bindings related to the exports of other packages, when your own custom package is loaded? If so, then why, except for convenience? It's just my opinion of course but in such a set up I would consider the behaviour of the global environment to be a strange thing, with behaviour that would be unusual and apparently inconsistent. It might seem useful in simple cases, but I would be concerned with addling the binding stack in such a way.

Can't you explicitly issue with(RealDomain): with(mypackage): at the top level instead, which keeps it clear what is loaded. Or, alternatively, could you create explicit exports for the names you want rebound by loading your package, even if those consisted simply of calls to exports of other packages such as RealDomain?

I do not see the behaviour in your first case (with a,b,c,d shown still) when running either Maple 16.00 32bit or Maple 16.02 32bit on ubuntu 10.04. I tried also on an i386 based SuSE 10.1 with 32bit Maple 16.00 and 16.02 and there again I also do not see the first behaviour you showed.

A Matrix (capital M) does not have last-name-evaluation. Even if you had used matrix (lowercase m) to get a table based array (which does have last-name-eval) then you should see the literal `A` printed as the result of issuing `A;` and not what you showed.

What minor version number are you using (Maple 16.00, 16.01, 16.02, etc)? Which Linux distribution? Could you double check whether you really do get that behaviour after issuing exactly the command shown?

Since you used uppercase M for Matrix, this is likely about evaluation under `print/rtable` I would expect, and not anything to do with last-name-evaluation.

acer

That call to Units:-AddSystem might be removed from the local AtomicWeight procedure. If you issued it just once at the top level, then the AtomicWeight procedure could be called repeatedly without its reexecuting that call to AddSystem each time. That should be more efficient. (I got a cpu timing difference of about 2 sec vs 24 sec, for 10000 repeats of a call like mat:-AtomicWeight(Cl).)

So then you could consider putting that call to Units:-AddSystem into a local module member named ModuleLoad. Once you save such a revised module to Library archive (the subject of your post) then that modification could come into play. In a new session (with the saved archive in libname) then that ModuleLoad procedure would get run just once at the moment that the module is first accessed from the archive. Ie, when you first call mat:-AtomicWeight or some other export, or when you first issue with(mat). It would also get called once upon the first reference following a restart.

Just a thought.

Apart from the section on modules in the Programming Manual, you might also take a quick look at this answer to a related, older question.

acer

@Gauss It's often best to state up front what you are really after.

If you know that the expression is just a product powers of names then you could be more direct, eg.

ee:=Pi^3*x^2*y^3/(z^7*r^3*p):

map(`[]`@op,select(type,[op(ee)],And(name,Non(constant))^negint));

                  [[z, -7], [r, -3], [p, -1]]

map(`[]`@op,[indets(ee,And(name,Non(constant))^negint)[]]);

                  [[p, -1], [r, -3], [z, -7]]

map2(op,1,[indets(ee,And(name,Non(constant))^negint)[]]);

                           [p, r, z]

And if you know that there will be no coefficients containing named constants then it can be simpler. Eg,

ee:=x^2*y^3/(z^7*r^3*p):

map(`[]`@op,select(type,[op(ee)],name^negint));

                  [[z, -7], [r, -3], [p, -1]]

map(`[]`@op,[indets(ee,name^negint)[]]);

                  [[p, -1], [r, -3], [z, -7]]

map2(op,1,[indets(ee,name^negint)[]]);

                           [p, r, z]

@cesar torres What happens if you change,

  A := Array(1 .. 5, 1 .. 3);

to

  A := Matrix(5, 3, datatype=float[8]);

?

What version of Maple are you using?

What is the value of `kk`?

acer

@Carl Love My mistakes, sorry. Indeed it is a very nice group effort. I have edited the attribution.

@taro yamada Which method you prefer may be a matter of taste, or depend on other as yet unmentioned restrictions.

One weakness of the first approach below is that it presumes knowledge about the equation: which side is a sum of terms, and which not. So it's not robust for blind re-use on other examples, say where the lhs and rhs are interchanged.

restart:
g:=x->int(x,k=0..M):
equ:=p(k)=alpha-gamma*q(k):
g(lhs(equ))=map(g,rhs(equ));

            int(p(k), k = 0 .. M) = alpha M + (int(-gamma q(k), k = 0 .. M))

restart:
g:=x->int(x,k=0..M):
f:=x->`if`(x::`+`,map(g,x),g(x)):
equ:=p(k)=alpha-gamma*q(k):
map(f,equ);

            int(p(k), k = 0 .. M) = alpha M + (int(-gamma q(k), k = 0 .. M))

restart:
define(g,linear,g(x::anything)='int'(x,k=0..M));
equ:=p(k)=alpha-gamma*q(k):
map(g,equ);

             int(p(k), k = 0 .. M) = alpha M - gamma (int(q(k), k = 0 .. M))

It's often better to state up front the specifics of your actual problem.

@taro yamada My earlier `map2` example did match your original description, because you wrote originally that you wanted map(f, lhs(x2)) and so your g:=x->int(x,k=0..M) would get mapped over the operands of the function call p(k). So it now seems that you only want `f` mapped across sums. Is that right?

I mean, I'm not sure that I understand what you want done to the right hand side of `equ`.

Do you want something other than the following?

equ:=p(k)=alpha-gamma*q(k);

                      p(k) = alpha - gamma q(k)

map(int,equ,k=0..M);

      int(p(k), k = 0 .. M) = int(alpha - gamma q(k), k = 0 .. M)

Or do you want the integral of the rhs sum to become a sum of integrals?

@taro yamada The Maple specific term is operands, and the mapping is done over the operands of a set or of an equation.

op( {a, b, c} );
                            a, b, c

op( y = a + b );
                            y, a + b

The `map2` command provides a convenience, and allows a shorter syntax than, say,

map( dummy->map(f,dummy), y=a+b );

                       f(y) = f(a) + f(b)

What has happened here (and in that `map2` example), is that `map` itself is mapped over the operands of the equation.

@Markiyan Hirnyk The question of low a degree a polynomial may approximate (to within a specific tolerance) the given expression on 0..1 is interesting theoretically, because of some potential numerical difficulties at both end-points.

In practice so-called range reduction (which is really domain reduction) could be used. And this allows the target error of 0.01 to be met by using only the restricted domain or 1/2..1 while computing the polynomial, since the mathematical function has a symmetry of reflection about the point (1/2,Pi/4) even if it lacks symmetry by reflection across a line.

If we take p as,

p:=1.135971747717227*x+.1950071341939759+.4815612713662410*(4*x-3)^3
-2.439038982375010*(4*x-3)^4-12.66183100870810*(4*x-3)^5
+41.62555613174249*(4*x-3)^6+164.1652491121847*(4*x-3)^7
+975.5300610456691*(4*x-3)^22+1785.055305607627*(4*x-3)^23
+11358.38836456872*(4*x-3)^18+23013.15164101551*(4*x-3)^19
+26106.22796515134*(4*x-3)^15-14347.77165088372*(4*x-3)^16
-30967.64393191083*(4*x-3)^17+.1043967698939514*(4*x-3)^2
-359.5777009952230*(4*x-3)^8-1176.044314414133*(4*x-3)^9
+1804.767818275263*(4*x-3)^10+5127.986643910732*(4*x-3)^11
-5625.866345976549*(4*x-3)^12-14311.88358476918*(4*x-3)^13
+11225.66738216834*(4*x-3)^14-5070.302621042330*(4*x-3)^20
-9728.730684483163*(4*x-3)^21:

then I believe that arcsin(sqrt(x)) may be approximated to within 0.01 over the doman 0..1 by using the translation evalf((Pi/2)-eval(p,x=1-x),16) for x in 0..1/2.

With Maple 17,

restart:
with(numapprox):
e:=proc(x) Digits:=100; arcsin(sqrt(x)); end proc:
c:=chebpade(e,1/2..1,[23,0]):
Digits:=16:
p:=sort(eval(c(x),T=orthopoly[T]));
sol:=piecewise(x>=1/2,p,evalf(Pi/2)-eval(p,x=1-x)):

Finding the minimal degree polynomial in the Chebyshev basis that meets the target tolerance for even the restricted domain seems tricky if using `numapprox`, and that might be due in part to some hard-coded constants within the `evalf/int/ccquad` procedure. I wonder whether some integrations could be done with more manual control.

@Preben Alsholm Does this get an error better than 0.016? Is it about 0.01543?

.785398163397445*T(0, 2*x-1)+.638238282327978*T(1, 2*x-1)
+0.723739716042111e-1*T(3, 2*x-1)+0.271442813438863e-1*T(5, 2*x-1)
+0.147364292749009e-1*T(7, 2*x-1)+0.969627362219072e-2*T(9, 2*x-1)
+0.722513270043944e-2*T(11, 2*x-1)+0.590232940948517e-2*T(13, 2*x-1)
+0.519621737331440e-2*T(15, 2*x-1):

Or how about,

.78539816339744830960*T(0, 2*x-1)+.63823828232797760017*T(1, 2*x-1)
+0.72373971604211358200e-1*T(3, 2*x-1)+0.27144281343886350513e-1*T(5, 2*x-1)
+0.14736429274900909502e-1*T(7, 2*x-1)+0.96962736221907198717e-2*T(9, 2*x-1)
+0.72251327004394449383e-2*T(11, 2*x-1)+0.59023294094851786867e-2*T(13, 2*x-1):

which might do better than 0.015?