PeridentityMatrix := proc(n::posint)
LinearAlgebra[LinearAlgebra:-CreatePermutation](Vector([
seq(n - i, i = 0 .. floor(1/2*n) - 1),
seq(floor(1/2*n) + 1 + i, i = 0 .. ceil(1/2*n) - 1)]));
end proc:
acer

fsolve should order the roots in a reproducible, non-session-dependent manner. This order can also be "fixed up" to agree with the ordering of indexed RootOf's.
But you might also want to make sure that your code is efficient. If you are trying to compute the floating-point approximations of several of the roots of a univariate polynomial, then one call to fsolve should suffice. On the other hand, evalf(RootOf(...,index=i)) could end up calling fsolve for each i. And that might be wasted, duplicate effort.
# Using your e1 and e2 as posted...
P:=convert(subs(E1=e1,E2=e2,E1*NZ^4-NZ^2*(2*E1^2)+E1^3-E1*E2^2=0),rational):
[seq(evalf(RootOf(P,NZ,index=i)),i=1..4)];
`RootOf/sort`([fsolve(P,NZ)]);
If you set stopat(fsolve) you can see that the first of these two approaches calls fsolve four times (each time computing all the roots) while the second approach calls fsolve just the once and then fixes up the ordering to agree with that of the indexed RootOf's.
acer

fsolve should order the roots in a reproducible, non-session-dependent manner. This order can also be "fixed up" to agree with the ordering of indexed RootOf's.
But you might also want to make sure that your code is efficient. If you are trying to compute the floating-point approximations of several of the roots of a univariate polynomial, then one call to fsolve should suffice. On the other hand, evalf(RootOf(...,index=i)) could end up calling fsolve for each i. And that might be wasted, duplicate effort.
# Using your e1 and e2 as posted...
P:=convert(subs(E1=e1,E2=e2,E1*NZ^4-NZ^2*(2*E1^2)+E1^3-E1*E2^2=0),rational):
[seq(evalf(RootOf(P,NZ,index=i)),i=1..4)];
`RootOf/sort`([fsolve(P,NZ)]);
If you set stopat(fsolve) you can see that the first of these two approaches calls fsolve four times (each time computing all the roots) while the second approach calls fsolve just the once and then fixes up the ordering to agree with that of the indexed RootOf's.
acer

It could be that the degree of estimated roundoff error is dependent upon the value of vm. As vm gets below some critical value, that estimation of roundoff might tip over some internal limit. After all, the integrand itself, and not just the upper limit of integration, depends on vm.
I doubt that evalf/Int would attempt to change epsilon from its initial value, in case of failure. More interesting is the question of whether it does, or should, change the working precision on the fly. I would say that it ought not to do so, aside from some fixed initial adjustment of addition of guard digits. It would be un-Maple-like. The degree of conditioning, or such, of a given problem would affect how high Digits might have to be to attain a given accuracy. In that case, there'd be some examples that could make evalf/Int go away for arbitrarily long periods of time, which few people might like.
But evalf/Int allows both epsilon and the working precision to be supplied as options. So there is control. It does what it's told to do. As mentioned, one can raise the working precision, hold epsilon constant, and attain solutions, with no method specified.
I found this interesting:
> foo := proc(v) if sqrt(sech(v) - sech(VM)) = 0 then return Float(undefined) \
> else return Re(cosh(v)/sqrt(sech(v) - sech(VM))) end if; end proc:
> N:=(vm,eps)->evalf(subs({VM=vm,EPS=ep\
> s},Int(eval(foo),0..VM,epsilon=EPS,method = _d01ajc))):
> infolevel[`evalf`]:=100:
> N(.000222,1.0e-6);
evalf/int/control: integrating on 0 .. .222e-3 the integrand
proc(v)
if sqrt(sech(v) - sech(0.000222)) = 0 then return Float(undefined)
else return Re(cosh(v)/sqrt(sech(v) - sech(0.000222)))
end if
end proc
evalf/int/control: tolerance = .10e-5; method = _d01ajc
Control: Entering NAGInt
trying d01ajc (nag_1d_quad_gen)
d01ajc: epsabs=.100000000000000e-8; epsrel=.10e-5; max_num_subint=200
d01ajc: procedure for evaluation is:
proc (v) if sqrt(sech(v)-sech(.222e-3)) = 0 then return Float(undefined) else \
return Re(cosh(v)/sqrt(sech(v)-sech(.222e-3))) end if end proc
Error, (in evalf/int) NE_QUAD_MAX_SUBDIV:
The maximum number of subdivisions has
been reached: max_num_subint = 200
Here, failure seems related to the maximal number of subdivisions allowed in the d01ajc integrator. But that doesn't seem to be a limit that is under user's control as an option to evalf/Int.
acer

One could add this: It is not at all uncommon for a function to have a singularity (or even something that merely appears as such) and be both very easy to symbolically integrate yet not so straightforward to numerically integrate. Not every numerical integration method will handle such functions the same way. So in evalf/Int, it's easily conceivable that method _Dexp could handle such an integrand without accumulating the same sort of accumulated roundoff error that method _d01ajc might accrue.
Without forcing a given method, evalf/Int can do analysis of the integrand and select whichever methods it wants to use to obtain the desired accuracy that is specified by the epsilon option.
acer

The integrand appears to have a singularity. Eg,
foo := proc(v) if sqrt(sech(v) - sech(VM)) = 0 then return Float(undefined) else return Re(cosh(v)/sqrt(sech(v) - sech(VM))) end if; end proc:
plot(subs({VM=1},eval(foo)),0.2..1);
So, if one attempts to integrate it numerically with method=_d01ajc, then many functional evaluations may be done near the singularity. That could account for enough accumulated roundoff error to prevent attaining an accuracy better than the smallest found successful epsilon.
For the plot, after,
N:=(vm,eps)->evalf(subs({VM=vm,EPS=eps},Int(eval(foo),0..VM,epsilon=EPS,method = _d01ajc))):
did you perhaps intend either of,
plot('N'(v,1.0e-6),v=0.2..1);
or,
plot(v->N(v,1.0e-6),0.2..1);
so as to delay any attempt at immediate evaluation of the argument N(v,1.0e-6) ?
acer

At a given working precision (Digits), the integrand might evaluate to a quantity with a possibly small nonzero imaginary component, for some value of v. That seems to be disallowed in the purely real integrator, Nag's d01ajc. So one could put a Re() call around the integrand, to help with the case of "unable to obtain a real result".
The other error message you obtained says that round-off error was apparently detected while trying to obtain the default tolerance. So one could try adjusting that, using evalf/Int's epsilon option.
Let's see how small one could make the tolerance, and have it succeed.
> foo := proc(v) if sqrt(sech(v) - sech(VM)) = 0 then return Float(undefined) else return Re(cosh(v)/sqrt(sech(v) - sech(VM))) end if; end proc:
> N:=(vm,eps)->evalf(subs({VM=vm,EPS=eps},Int(eval(foo),0..VM,epsilon=EPS,method = _d01ajc))):
> N(.01,1.0e-6);
2.221566392
> N(.01,1.0e-7);
2.221566471
> N(.01,1.0e-8);
2.221566424
> N(.01,1.0e-9);
Error, (in evalf/int) NE_QUAD_ROUNDOFF_EXTRAPL:
Round-off error is detected during extrapolation.
> Digits := 15:
> N(.01,1.0e-9);
Error, (in evalf/int) NE_QUAD_ROUNDOFF_EXTRAPL:
Round-off error is detected during extrapolation.
I didn't expect setting Digits to 14 or 15 to make much, if any, difference, as that's going to be a close match to a compiled double precision external integrator using evalhf callbacks for functional evaluations of the integrand.
> Digits := 25:
> N(.01,1.0e-9);
Error, (in evalf/int) expect Digits<=15 for method _d01ajc but received 28
This error message seem to be indicating that it's not valid to use this method with a working precision greater than evalhf's.
I notice that if one removes the method=_d01ajc restriction then at Digits=25 a result is returned even with epsilon as small as 1.0e-9 .
I do not understand the remark, "Curiously, these problems occur where elementary integration is feasible." Numerical computational issues can have nothing to do whether a problem is symbolically solvable, in general and for quadrature in particular, as an entirely different methodology might be used. I don't see what is remarkable about such a difference.
acer

One way to experiment with whether the extra time may be due to garbage collection (gc), is to increase the limit that controls gc frequency.
For example, set,
> kernelopts(gcfreq=10^8):
and then see what efect that has on the repeated loop timings.
acer

The original example is covered by its own explanatory paragraph in the help-page ?sign .
There is also sometimes confusion between the functions sign() and signum().
Consider this below, where no quoting is necessary,
plot([signum(x)],x=-1..1);
and see ?signum .
acer

Hello again,
The Matrix ls_m is currently set up as 6x9. So, assuming that the extra 3 columns of ls_m actually contain pertinent and meaningful data, the system appears underdetermined. Presumably, this is why you've supplied a name for the generation of the free variables to appear in any solution, by supplying the option free='t'.
What might not be made clear enough in the help-pages ?LinearSolve and ?IterativeSolver is that the underlying sparse solvers are not designed for producing symbolically parametrized solutions of underdetermined systems. They are purely numerical, and can by nature only ever generate a purely floating-point solution with no free variables. (Such iterative methods generally use some form of Matrix-Vector multiplication, or a related way to "iteratively" get closer to an answer. There's no room in that for supplying any symbolic piece to the result.)
If ls_m were 6x6 (or 9x9) and had datatype=float[8], with B of corresponding size, then that should work with the sparse method, I believe. You might even be able to pad out the data so as to get the a (single, purely numeric) sol which approximately satisfied the original system. But the current sparse solvers are purely numeric, and can't get you a parametrized solution, I don't think.
As you've pointed out, though, this below does produce a solution to the problem (sized, as stated), as long as ls_m and B don't contain floats,
sol := LinearSolve(ls_m, B, free='t');
So perhaps you could just use that, then, and conclude that the sparse solver is float-based and not appropriate here?
acer

Hello again,
The Matrix ls_m is currently set up as 6x9. So, assuming that the extra 3 columns of ls_m actually contain pertinent and meaningful data, the system appears underdetermined. Presumably, this is why you've supplied a name for the generation of the free variables to appear in any solution, by supplying the option free='t'.
What might not be made clear enough in the help-pages ?LinearSolve and ?IterativeSolver is that the underlying sparse solvers are not designed for producing symbolically parametrized solutions of underdetermined systems. They are purely numerical, and can by nature only ever generate a purely floating-point solution with no free variables. (Such iterative methods generally use some form of Matrix-Vector multiplication, or a related way to "iteratively" get closer to an answer. There's no room in that for supplying any symbolic piece to the result.)
If ls_m were 6x6 (or 9x9) and had datatype=float[8], with B of corresponding size, then that should work with the sparse method, I believe. You might even be able to pad out the data so as to get the a (single, purely numeric) sol which approximately satisfied the original system. But the current sparse solvers are purely numeric, and can't get you a parametrized solution, I don't think.
As you've pointed out, though, this below does produce a solution to the problem (sized, as stated), as long as ls_m and B don't contain floats,
sol := LinearSolve(ls_m, B, free='t');
So perhaps you could just use that, then, and conclude that the sparse solver is float-based and not appropriate here?
acer

You wrote of transmission0 and reflection0, that, "all they contain are simple floats."
As described, however, the code is generating a great many so-called software floats. Maple will dispose of those, once no longer needed, during garbage collection. But that extra garbage collection adds to the total overhead. If you really just need the Matrices to store floating-point numbers, and if hardware double-precision is sufficient, then consider placing an appropriate datatype on the Matrices.
For example,
> transmission0:=Matrix(1200,2,datatype=float[8]):
and similarly for reflection0, re, and te.
You might also be able to fill the first column of both transmission0 and reflection0 more quickly, by pulling the assignments (into their first columns) out of the loop. Completely outside of the loop, create a Vector,
> firstcol := Vector(1200,(j)->0.0416666666*j,datatype=float[8]):
and also do, outside the loop,
> transmission0[1..1200,1] := firstcol:
> reflection0[1..1200,1] := firstcol:
I hope I got that right.
There might be more efficient ways to initialize the second column too, but it's not possible to say without knowing how rcwa writes to te and re depending on M and how M changes with m.
Also, if you made rcwa accept re and te as parameters (it should still be able to update them efficiently, inplace, if they had float[8] datatype), while also making rcwa return M instead of writing to it as a global, then perhaps you could then run rcwa under the faster purely floating-point evalhf interpreter. See the help-page ?evalhf for details.
acer

I am curious why there is interest shown here in speed performance but not in memory performance. For example, how do the various techniques differ in terms of how much extra memory they allocate? I find that Maple performance discussions tend in general to have an emphasis about speed of execution. Isn't a performance metric which combines both speed and memory allocation just as valuable, and perhaps more valuable in general?
I have seen a few different performance metrics commonly used in the past. One is time_duration*allocation_increase, and another is time_duration*(allocation_increase^2).
I prefer the technique Roman used here, to have restarts between comparisons of various methods. For, even if one is not measuring memory allocation as part of the performance metric, there can be a penalty for garbage collection which changes with ongoing differences in memory allocation. So it can happen that the order of tests affects how they perform. Putting distinct tests into distinct Maple runs can alleviate such obscuring interference. (As we know, it was not always true as it is in Maple 10 that memory would be completely freed upon restart. So multiple sessions for distinct tests is cleaner all round.)
acer

Despite the fact that eq contains coefficients whose exponents grow (exponent of 10) with respect to evaluation at Digits precision, the numerical solution appears to remain stably computed. So perhaps the solution is accurate.
The left-hand-side of eq appears to be monotonic and continuous from -10^8 to 10^3. (Wouldn't it be mildly surprising if the *exact* solution of equations which well-represent beam mechanics not at least be continuous?) You mentioned that you had a class of equations to handle. Will they all be monotonic, or at least continuous? If so, maybe it could suffice to find a value of N which makes lhs(eq) positive and another which makes it negative, in which case you know that there must exist a value between them which satisfies eq. You could try peppering a numeric range for N, in an attempt to find such a pair of values of N. That'd save the effort of asking fsolve to locate and closely approximate an actual root.
As far as the code goes, the final eq and solution appear to remain unchanged if you do the following:
- replace the subs() calls with eval() calls
- replace the evalf(...int()) calls with evalf(...Int()) calls
- replace the int() call when assigning to f by an Int() call
This seems to speed it up.
acer

Despite the fact that eq contains coefficients whose exponents grow (exponent of 10) with respect to evaluation at Digits precision, the numerical solution appears to remain stably computed. So perhaps the solution is accurate.
The left-hand-side of eq appears to be monotonic and continuous from -10^8 to 10^3. (Wouldn't it be mildly surprising if the *exact* solution of equations which well-represent beam mechanics not at least be continuous?) You mentioned that you had a class of equations to handle. Will they all be monotonic, or at least continuous? If so, maybe it could suffice to find a value of N which makes lhs(eq) positive and another which makes it negative, in which case you know that there must exist a value between them which satisfies eq. You could try peppering a numeric range for N, in an attempt to find such a pair of values of N. That'd save the effort of asking fsolve to locate and closely approximate an actual root.
As far as the code goes, the final eq and solution appear to remain unchanged if you do the following:
- replace the subs() calls with eval() calls
- replace the evalf(...int()) calls with evalf(...Int()) calls
- replace the int() call when assigning to f by an Int() call
This seems to speed it up.
acer