@Markiyan Hirnyk I have not so far stated that there are only two pairs of purely real solutions.

But I think that you should at least look harder at that 3rd result you cited from DirectSearch, where `a` was approx. 15.2 and `x` was approx. 5e-7.

For example, why would eq[2] become less than 2 when `x` gets close to zero?

Also, if one creates two nonevalhable procedures to evaluate eq[1] and eq[2] while respecting a very high working precision then I don't see the solution curves crossing near x=0.

For example, zooming in...

U1:=unapply(eq[1],[a,x]):
F1:=proc(a,x) Digits:=500; []; U1(a,x); end proc:
U2:=unapply(eq[2],[a,x]):
F2:=proc(a,x) Digits:=500; []; U2(a,x); end proc:

plots:-implicitplot( [F1,F2], 33..35, 0..1e-14,
                     labels=[a,x], color=[red,blue], gridrefine=2 );

In other words, I don't believe that the red and blue curves are going to cross at some larger value for `a`, for the curves seen going off to the bottom right of this,

plots:-implicitplot( [F1,F2], 1..5, 0..0.5,
                     labels=[a,x], color=[red,blue], gridrefine=2 );

If they do seem to cross, then zooming in (and also possibly increasing precision or grid refinement) seems to show such a manifested crossing as spurious. Of course some deft exact analysis (or the difference eq[2]-eq[1], perhaps?) would be better.

@awass The key difference is not that your example was (or was not) inside a do-loop. The key difference was that one of your examples contained literal numerics which automatically simplified to zero in a denominator. (In your looped version one of the literal 7s was replaced by the loop index j.)

Automatic simplification precedes usual evaluation, and as a consequence of which is not delayed by unevaluation quotes and not trapped by the try..catch mechanism. In other words you cannot disable it for your particular example except by using something other than literal numeric input.

By the way, some numeric exception in your examples can also be managed by another mechanism -- a custom Numeric Event Handler. Eg,

restart:

'1/(7-7)';
Error, numeric exception: division by zero

NumericStatus( ':-division_by_zero' );
                              true

MyHandler := proc(operator,operands,default_value)
   NumericStatus( division_by_zero = false );
   return infinity;
end proc:

NumericEventHandler( ':-division_by_zero'=MyHandler ):

1/(7-7);
                            infinity

NumericStatus( ':-division_by_zero' );
                             false

That's just an example of usage. You could alter it to return a different value, or alter it to not reset the numeric event's status flag, etc.

There are other kinds of automatic simplification that can occur than just the purely numeric subtraction & division example you gave. Eg,

'2*(x+y)';

                             2 x + 2 y

'proc() x*~y; end proc';

                 proc() ~[:-`*`](x,  $, y) end proc;

'x-x';

                                 0

@Stephen Forrest Compiling a symbolic piecewise produced by CurveFitting:-Spline is not a good approach in the case that there are many data points. The task of producing the piecewise scales badly with the number of points, and the compilation cost must be paid for each new set of data. That compiled-piecewise method is easy to set up, but its set up is very costly at large sizes. This is the kind of scenario for which ArrayInterpolation can be used for much more efficient set up. The addition of an even easier way to set up a proc that does lean scalar->scalar calls to ArrayInterpolation would improve Maple.

@Andriy I would concur with what Preben has said about `A` not being evaluated when say `A[1]` is called, and that is key.

If you are going to access elements of `A` many times then you would be better off assigning Q:=eval(A) after the assignment b:=a, and then working with `Q`. 

Note that if you had done all this inside a procedure with local `A` then numelems(A) would return 3 even after the assignment of `a` to `b`. That would be similar to calling numelems(eval(A,1)) at the top level which Preben described. In that procedure setting the extra statement like Q:=eval(A) might seem less unusual.

It might seem that another reasonable workaround for the top level is to use op(2,A) instead of A[2]. That should produce the results you were expecting, ie. op(1,A) returns `a` and op(2,A) returns `c`. But I suspect that for a large set `A` with many elements there could be large overhead with using `op` instead of square bracket indexing to handle this special situation that some entries change upon evaluation. I would worry that each new expression sequence of evaluated entries might get resorted during the uniquification process, upon each (higher than 1-level) evaluation of A.

@Markiyan Hirnyk I used the example about which the original poster had inquired. You're not contributing in a meaningful way here.

@sbh 

It worked for me in Maple 17.

restart:
main2 := X^Q = V^P*(V-1)^Q/(c^(P-Q)*(-c^2+V)^Q):
split := X^Q = V^P*(a[1]+a[2]*V+a[3]*V^2+O(V^3)):
one := (V-1)^Q/(c^(P-Q)*(-c^2+V)^Q):
N := 5:
two := simplify(series(one,V,N)) assuming real, V>c^2, c>1:
three := convert(two,polynom):
four := simplify(V^P*three,power):
ans1 := frontend(collect,[four,[seq(V^(P+i),i=1..N)]],[{`*`,`+`,list},{}]);
frontend(collect,[four,[seq(V^(P+i),i=1..N)],
                        u->factor(simplify(expand(u*c^(P+Q)))
                                           /c^(P+Q))],
         [{`*`,`+`,list},{}]):
ans2 := sort(eval(%));
simplify(ans1 - four);
simplify(ans1 - ans2);

I think that such a point plot which displays the ordered complex roots of a univariate polynomial (alongside each index value) would be a useful addition to the ?RootOf/indexed help-page. By that I mean useful for students.

@Markiyan Hirnyk Your plot misses the aspect of which red dot corresponds to which indexed RootOf, which to the new user may provide more insight.

I was trying to illustrate not just where the index=3 RootOf would appear (on the real line) but also how the complex roots get sorted altogether. Hence I made the extra effort in the plot and printing, going as far as writing extra code to ensure that a new user might get the whole concept more clearly (in particular, from the visual plot alone).

It's a true shame that this site (or its maplenet backend?) ruins inlined worksheets by automatically adding spurious gridlines on 2D plots.

@Thomas Richard Upon reading the Post it occured to me that perhaps the .bin installer file might not be set as executable.

@Markiyan Hirnyk Please tell us the point of your remark, because without any words of explanation it seems pointless. What were you hoping that your command would return, in terms of an explicit, exact symbolic result?

@Andriy Do you mean that you now want to unassign the names in `a`?

restart:

a:=[chris,jack,john,jason,fred,alex]:
seq(assign(i, 0), i = a);

map(unassign, eval(a,1))[];

a;
             [chris, jack, john, jason, fred, alex]

@J4James I only put the plottools:-transform in there to reflect the x- and y-axes. That command turns the data from Array into listlist. So either use some invocation involving nops (or numelems) instead of rtable_dims to pick off the dimensions of the two listlists, or don't transform the axes at that time.You don't necessarily need to transform and switch x-y axes, as it is easy to just change how the entries get printed.

F12 := -(1/2*(S+sqrt(S^2+4*alpha)))*alpha:
F22 := (1/2*(-S+sqrt(S^2+4*alpha)))*alpha:
p1 := plot3d({F12,F22}, alpha=max(-1,-S^2/4)..0, S=-20..20):
A:=op([1,1],p1):
m:= rhs(rtable_dims(A)[1]):
n:= rhs(rtable_dims(A)[2]):
for i from 1 to m do
for j from 1 to n do
X := A[i,j,1];
Y := A[i,j,2];
Z := A[i,j,3];
printf("%f %f %f\n",Y,X,Z);
end do:
end do:
B:=op([2,1],p1):
m:= rhs(rtable_dims(B)[1]):
n:= rhs(rtable_dims(B)[2]):
for i from 1 to m do
for j from 1 to n do
X := B[i,j,1];
Y := B[i,j,2];
Z := B[i,j,3];
printf("%f %f %f\n",Y,X,Z);
end do:
end do:

@J4James Yes.

op(p1) are all the operands of the PLOT3D structure assigned to `p1`. It has two MESH substructures, because the first argument to plot3d was a set of two expressions (lambda1 and lambda2).

op(1,p1) is the first MESH, and op(2,p1) is the second MESH.

And so op([1,1],p1) is the Array in the first MESH, and op([2,1],p1) is the Array in the second MESH.

@J4James It appears that you have some preconception of how the data that may represent a surface must be formatted or structured. There is no single way in which the data that represents a surface (to somebody or to some application) must be structured.

Are you saying that you need or expect a regular grid of data points -- evenly spaced in x and y directions say -- because you wish to use the data inside some other (plotting?) program which expects such a format? If so, then you could make that more clear, if that is what you mean by "extracting" the data.

What exactly do you consider to be wrong with the extracted data (triples) from your original posted code? Is it because some other application also renders one curved edge as very jagged? If that is so, then what format would you want instead? Using a very, very fine grid just to get rid of jaggedness that could otherise be far more easily handled by a different structure&interpretation seems like a resource expensive approach.

What is your target application for the exported data? What is the defintion of the format in which you want data to be exported?

The plottools:-getdata command is just an alternative (to `op`) for getting one's hands on the Arrays or data. It doesn't provide anything different about how that data is formatted. (All it really does is call `op` itself).

lambda1:=(S+sqrt(S^2+4*alpha))/(2):
lambda2:=(S-sqrt(S^2+4*alpha))/(2):
p1:= plot3d({lambda1,lambda2}, alpha=max(-1/4*S^2,-5)..5, S=-10..10, grid=[5,5]):
A:=op([1,1],p1):
B:=op([2,1],p1):
dd:=plottools:-getdata(p1);
dd[1,-1]; # A
dd[2,-1]; # B

The joy of the MESH structure for your example is that Maple's interfaces know how to render a surface from it that doesn't have the very jagged edge. There's no reason to expect that some other application should be able to do the same with it.

@J4James 

lambda1:=(S+sqrt(S^2+4*alpha))/(2):
lambda2:=(S-sqrt(S^2+4*alpha))/(2):
p1:= plot3d({lambda1,lambda2}, alpha=max(-1/4*S^2,-5)..5, S=-10..10, grid=[5,5]):
A:=op([1,1],p1):
m:= rhs(rtable_dims(A)[1]):
n:= rhs(rtable_dims(A)[2]):
for i from 1 to m do
for j from 1 to n do
X := A[i,j,1];
Y := A[i,j,2];
Z := A[i,j,3];
printf("%f %f %f\n",Y,X,Z);
end do:
end do:

Make sure that `A` is assigned a MESH Array. And repeat as necessary, for other MESH Arrays in the PLOT3D structure.

In modern Maple plottools:-getdata can be used instead of `op`, for getting at the rtables or listlists in the MESH or GRID. But that still leaves you the work of exporting it.