@Mac Dude It seems like a bug in ScientificConstants, where it doesn't properly accomodate the new system.

This next looks ok,

restart:

with(ScientificConstants):
em:=Constant(electron_mass):

GetValue(em),GetUnit(em);

                            -31                   
              9.109381882 10   , Units:-Unit('kg')

But now, with a system with energy and action to be simplified in terms of MeV and MeV*s respectively,

restart:

Units:-AddSystem('Accelerator',Units:-GetSystem('SI'),MeV,MeV*s); 
Units:-UseSystem('Accelerator');

with(ScientificConstants):
em:=Constant(electron_mass):

GetValue(em),GetUnit(em);

                            -18                   
              5.685626500 10   , Units:-Unit('kg')

@Mac Dude It seems like a bug in ScientificConstants, where it doesn't properly accomodate the new system.

This next looks ok,

restart:

with(ScientificConstants):
em:=Constant(electron_mass):

GetValue(em),GetUnit(em);

                            -31                   
              9.109381882 10   , Units:-Unit('kg')

But now, with a system with energy and action to be simplified in terms of MeV and MeV*s respectively,

restart:

Units:-AddSystem('Accelerator',Units:-GetSystem('SI'),MeV,MeV*s); 
Units:-UseSystem('Accelerator');

with(ScientificConstants):
em:=Constant(electron_mass):

GetValue(em),GetUnit(em);

                            -18                   
              5.685626500 10   , Units:-Unit('kg')

Answers about how to do this best (or perhaps just better) may depend on the particular nature of `a1`. Could you provide some details in the form of a fully functioning, explicit example?

acer

@Markiyan Hirnyk I interpreted the question as being about two things: the holes in the plot, and the long computation time.

Maple can get quirky in strange ways when Digits is set less than 5, so having it that low is not a great idea.

The holes in the plot are because the individual quadrature attempts failed (for input pairs in the plane). evalf/Int infers the value for epsilon based on Digits, when that option is supplied. But if Digits is set too low then the inferred looser tolerance may not be helpful as there might bot be enough working precision even to satisfy the looser tolerance. Hence quite often it helps to have a higher working precision (default Digits might do) while forcing a looser tolerance separately.

And evalf/Int might converge enough to satisfy a looser tolerance more quickly. Hence I suggested leaving Digits at its default value (10) while supplying looser tolerance (larger epsilon). That appears to help with both the failing values as well as the speed.

Having Digits be as high as 8 (or 10, the default) might fix the holes. But reducing Digits from default 10 down to a value of 8 doesn't do as good a job with speeding it all up as does supplying a looser tolerance. So typing in Digits:=8 or what have you just seems to be unnecessary typing, that might also obscure what may matter more. I did not test whether Digits=10 is needed; it might be, but even if not there seems less reason to reduce it from its default as doing so does not by itself cure the speed issues.

And the same may go for other tweaks such as using the non-iterated quadrature method and reducing the plot points.

@Markiyan Hirnyk I interpreted the question as being about two things: the holes in the plot, and the long computation time.

Maple can get quirky in strange ways when Digits is set less than 5, so having it that low is not a great idea.

The holes in the plot are because the individual quadrature attempts failed (for input pairs in the plane). evalf/Int infers the value for epsilon based on Digits, when that option is supplied. But if Digits is set too low then the inferred looser tolerance may not be helpful as there might bot be enough working precision even to satisfy the looser tolerance. Hence quite often it helps to have a higher working precision (default Digits might do) while forcing a looser tolerance separately.

And evalf/Int might converge enough to satisfy a looser tolerance more quickly. Hence I suggested leaving Digits at its default value (10) while supplying looser tolerance (larger epsilon). That appears to help with both the failing values as well as the speed.

Having Digits be as high as 8 (or 10, the default) might fix the holes. But reducing Digits from default 10 down to a value of 8 doesn't do as good a job with speeding it all up as does supplying a looser tolerance. So typing in Digits:=8 or what have you just seems to be unnecessary typing, that might also obscure what may matter more. I did not test whether Digits=10 is needed; it might be, but even if not there seems less reason to reduce it from its default as doing so does not by itself cure the speed issues.

And the same may go for other tweaks such as using the non-iterated quadrature method and reducing the plot points.

I googled kovacic algorithm and the first hit was this and most of the first page of hits seemed relevant. I mention that first hit because its references indicate a preprint (1979?) by Kovacic as well as an implementation by D.Saunders given at ACM 1981.

Another hit that stood out was this Maple help-page (even if it may differ in implementation), which cites Kovacic at end in its References section,

acer

@Carl Love I was not claiming that it is bug in 2-argument eval. I stated that `eval/if` relies on the behaviour of 2-argument eval.

I'm not sure that `eval/if` is quite right. There are other corners, too.

The routine `eval/if` is affected by the following behaviour of 2-argument `eval`,

> eval('sin(r)', r=0);

                                       0

> eval('sin(0)', r=0);

                                    sin(0)

Note that `seq` does not behave like that,

> seq('sin(r)', r=0); 

                                       0

> seq('sin(0)', r=0);

                                       0

acer

@Preben Alsholm Yes, thanks, that's why I included my second example. It's an oddity amongst oddities.

> restart:

> eval(`if`(r,sin(Pi),p),r=true); # hmm

                                    sin(Pi)

> f:=x->x:

> eval(`if`(r,f(r),p),r=true);

                                     true

> eval(`if`(r,f(2),p),r=true); # hmm

                                     f(2)

> seq(`if`(r,f(2),p),r=true);

                                       2

It looks like a bug.

acer

What is the purpose of the first loop in the SampleQ procedure (that loops for i from 1 to Size)?

acer

@AndreaAlp It wasn't clear to me what was your exact problem. I thought that perhaps you wanted emitted Matlab code that contained elementwise ~ operators, but weren't getting them.

Perhaps you could upload a short but representative example, so that it would be more clear what you need.

@AndreaAlp It wasn't clear to me what was your exact problem. I thought that perhaps you wanted emitted Matlab code that contained elementwise ~ operators, but weren't getting them.

Perhaps you could upload a short but representative example, so that it would be more clear what you need.

I don't think that it is clear what you mean by, "Since g(x,y) has stationary  point at (0,0), I have to interpolate f(x,y) except (0,0) but also quite close to (0,0)."

Could you please justify that? Why do you believe that numeric quadrature would fail? Why do you believe that only a symbolic interpolatory result would suffice? And what would you do with such a result?

If you only intend to use a final symbolic interpolant for plotting or other computational tasks then an assertion that only a symbolic interpolant would suffice seems suspect. A numeric black box interpolating function might also do.

It's not even clear why you must interpolate, if you can poll your `f` and `g` at any x-y pair.

acer

I may be missing something, but is there some significant way (performance, or other) in which this differs from the splitting by and divvying up with the Threads:-Task model, to parallelize? (See the attachment in my worksheet in the comment to your earlier response above, where I did that.) It seemed to me that both get a small bit of speedup, but not a by large factor.

...Almost as if some Maple level overhead was dealt with better (parallelized, perhaps) while a large significant portion of the computation (the externally called bit perhaps) was done in the same manner.

And CodeTools:-Usage shows a similar value for both the wallclock time[real] and the cpu time (summed for all threads), which indicates that the thing may not actually running concurrently all of the code done under the maple threads. See my comments about possible mututal blocking by call_externals -- which is not confirmed.

It is also possible that multiple cores (not maple threads) may already have been doing some of the external work in parallel, via BLAS.

acer