Simply construct, as the objective to be fed to the optimizing command, a procedure that takes a single argument like `t` and internally calls x(t) and y(t) before calling and returning the value of f(x(t),y(t),t).

obj := proc(t) f(x(t),y(t),t); end proc;

TheOptimizer( obj, ... );

Or are you trying to say that you want to pass `x` and `y` as well, in the call to TheOptimizer? One easy way is to use single-right quotes, which are usually called uneval quotes because they delay evaluation.

obj:=proc(T,X,Y) f(X(T),Y(T),T); end proc:

f:=proc(a,b,c) a*c + b; end proc:

x := z->z^2-10*z+1:

y := z->sqrt(z):

Optimization:-Minimize( 'obj'(t, eval(x), eval(y) ), t=2..10 );

                   [-138.93535451785112, [t = 6.606395251725914]]

#plot( 'obj'(t, x, y), t=2..10 );

acer

The paradigm of producing all frames of an animation "up front" is old and dusty, and ought to be retired from Maple.

It is now possible to produce the majority of both 2D and 3D plots so quickly and efficiently that a frame-on-demand paradigm of animations should be ushered in. The PlotComponent is the vehicle for getting this.

Part of the speed is due to machines (and thus Maple's kernel/Library runtime) becoming faster. Part is due to relatively recent developments such as being able to use (and populate) hardware float precision 'float[8]' rtables for the data portions of plot data structures.

In contrast, the GUI has become more and more bloated (and leaky!). Fifteen years ago one could often manage not much more than a 25-50 frame animation, on account of the time it took to produce all frames up front as well as the total memory cost required to store all frames at once. Amazingly, the situation is too often much the same today!

The idea that one can raise the old-style animation's frame-per-second play rate is a bit silly. If you can only bearably produce 25-100 frames then a 60-frame-per-second play rate is unimpressive. Blink twice and the movie is finished. And, as many people have experienced, watching the Standard GUI consume 700MB of memory to produce and play a 100-frame old-style animation is disheartening, but subsequently watching it allocate another 700MB and freeze up while exporting that same animation is futility.

In contrast, many frame-on-demand new-style "animations" can run to hundreds and hundreds of frames while keeping memory down below 50MB. The fastest techniques of populating the plot data structures means that a 30-40 frame-per-second play rate is often possible (and it can be slowed down, with an added variable-length delay).

The Application worksheet mentioned in the parent Post has four 50-frame old-style animations. It took a staggering 60 seconds (!!) and about 380MB allocated main memory for a fast i7 to load that worksheet and autoexecute code (which produces all the frames of the animations). That is... a long time.

I replaced all four animations in the worksheet with a triple of Plot, Slider, and Label Embedded Components. All the action goes inside the Slider. I put all four of those triples inside that Application workseet, and removed the old auto-executed animations (but retained the 2 procedures which make the plot frames). The new version loads immediately. It stays at 50MB memory allocated, on a 64bit machine. The sliding action is slower, but then I haven't bothered to improve its plotting frame-production code. Download Paratry.mw

This is by far the first time that these issues have come up on this site. But there is now little point in arguing about the basic premise any more. What is worth discussing is the need for accompanying support. The PlotComponent currently has about a dozen programmable properties related to old-style animations. Those are a waste of space. (I would give them all up for programmable theta/phi/psi orientation properties.) What is needed, moving forwards, is an new-style animation frame-capture and export facility from PlotComponents. Preferably something more modern than animated GIF. Even mpeg-2 would be better.

Don't misunderstand: there is still work to be done. Getting a fast enough frame production can require techniques such a re-using the underlying float[8] rtable (Joe Riel shows the effectiveness of that here). But it's worth doing. (Did you know that Mathematica got a mostly-hidden code-compiler way back in release 2.2 or something, and that it can get used for auto-compiling expressions to be plotted? Maple could do that too, and hook it up to a super fast float[8] plot-rtable populator.)

acer

Yes, there were some changes to `solve` between 10.00 and 11.00, more specifically that appeared in some of the point-releases (10.01 through 10.06), if I recall correctly.

Now, I don't remember the specifics offhand, but if I do recall that there were some issues (such as you describe) introduced with some of those changes. I also believe that at least some of those have been remedied in the past few years (ie, with releases 14 and 15).

Can you upload any worksheet with some example of a problematic system, so that at the very least we might test and check whether in fact such problems persist in the current (15) major release?

One thing you could try is to call a polynomial system solver via an entry point like SolveTools:-PolynomialSystem. That might help if some of the problem is in preliminary effort made under `solve`, before despatch to such a routine. (Note that the Maple 11 help-page for that routine is quite different from the more enhanced version in Maple 15, whose help-page is online here.) Yet another possibility might be to try and leverage calling the Groebner package's routines directly.

acer

Your objective is an equation, which makes no sense. One can optimize a scalar value.

Did you instead intend to use lhs(eqns)-rhs(eqn)?

acer

In modern Maple, producing a plot can boil down to populating a (datatype=float[8]) Array or Matrix with computed data. My most recent experiences with this have been that if that work can be done entirely within evalhf then yes, it can be multithreaded. Often quite easily, with the Task threading model.

Also, temporarily escaping evalhf by using `eval` to call-back to regular Maple doesn't cut it here; it has to be done entirely withing evalhf.

However, most useful things in the Maple Library are not thread-safe when run outside of evalhf mode. And lots of things cannot be run under evalhf. Examples include `dsolve`, `is`, and `int`, etc, etc.

Moreover, a great deal of what can be run under evalhf can all be Compiled. Some evidence indicates to me that a non-threaded (serial) Compiled implementation can edge out a quad-core threaded (parallel) evalhf implementation of populating a float[8] Array.

nb. This is why a thread-safe runtime of the Compiler is needed, as that would be all round, significantly fastest.

acer

Another approach might be to set it all up as a single coupled system, instead of iterating the dsolve/numeric process.

Relative performance seems to depend on N.  At N=10 it is twice as fast as iterating (and at N=8 thrice as fast), provided that you only want to compute values for the last y[i](t). When plotting all individual y[i](t) together, the better performance is not seen.

At higher N=16, both steps take long, and the benefit has vanished. Can it be improved? 

Re-using the example from elsewhere (but getting rid of the `abs`).

First, with the iterated approach,

restart:
st:=time():
sys:= diff(y(t),t) = cos( 1/2*y(t)-x(t)-Pi^2 ):
ini:= y(0)=10*sqrt(2)/2:
formula := (FF(t))^2 - sin(y(t))^2 - cos(y(t)) - 0.2*sqrt(t):
F[0]:=sin:
# This approach starts counting from 1
N:=8:
for i from 1 to N do
   newsys:=subs(x=F[i-1],sys);
   sol[i]:=dsolve([newsys, y(0)=10*sqrt(2)/2], numeric,
                  output=listprocedure, known=F[i-1]):
   eval(y(t),sol[i])(0.2);
   F[i]:=unapply('evalf'(subs(y(t)=subsop(3=remember,eval(y(t),sol[i]))(t),
                              FF=F[i-1],
                              formula)),
                 t,numeric,proc_options=[remember]):
end do:
solvetime := time()-st;
#st:=time():
#plot([seq(F[k],k=1..N)],0.1..2.0,
#     legend=[seq(F[k],k=1..N)],
#     color=[seq(RGB(0,1/((N-k+2)/7)^2,1-1/((N-k+2)/2)),k=1..N)],
#     thickness=2);
#time()-st;
st:=time():
plot([seq(F[k],k=N..N)],0.1..2.0,
     legend=[seq(F[k],k=N..N)],
     color=[seq(RGB(0,1/((N-k+2)/7)^2,1-1/((N-k+2)/2)),k=1..N)],
     thickness=2);
runtime:=time()-st:
solvetime + runtime;

And now as a coupled system.

restart:
st:=time():
sys:= diff(y[i](t),t) = cos( 1/2*y[i](t)-x(t)-Pi^2 ):
ini:= y(0)=10*sqrt(2)/2:
formula := (FF)^2 - sin(y(t))^2 - cos(y(t)) - 0.2*sqrt(t):
F[1]:=sin(t):
newsys[1]:=subs(x(t)=F[1],sys):
# This approach starts counting from 1
N:=9:
for i from 2 to N do
   F[i]:=subs(y(t)=y[i-1](t),FF=F[i-1],formula);
   newsys[i]:=subs(x(t)=F[i],sys);
end do:
sol:=dsolve([seq(newsys[i],i=1..N),
                seq(y[i](0)=10*sqrt(2)/2,i=1..N)], numeric,
              output=listprocedure):
allysols:=seq(y[i]=eval(y[i](t),sol),i=1..N):
allfuncs:=[seq(unapply(subs(allysols,F[k]),t,proc_options=[remember]),
              k=1..N)]:
seq(allfuncs[i](0.02),i=1..N):
solvetime := time()-st;
#st:=time():
#plot([seq(subs(allysols,F[k]),k=2..N)],
#     t=0.1..2.0,
#     legend=[seq(eval(y[k](t),sol),k=2..N)],
#     color=[seq(RGB(0,1/((N-k+2)/7)^2,1-1/((N-k+2)/2)),k=2..N)],
#     thickness=2);
#time()-st;
st:=time():
plot([seq(subs(allysols,F[k]),k=N..N)],
     t=0.1..2.0,
     legend=[seq(y[k](t),k=N..N)],
     color=[seq(RGB(0,1/((N-k+2)/7)^2,1-1/((N-k+2)/2)),k=1..N)],
     thickness=2);
runtime:=time()-st:
solvetime + runtime;

So, you might experiment with your larger system. The relative aspects of the performance might vary with example.

acer

M:=[(s-s1)*(s-s2), (s-s0)*(s-s2), (s-s0)*(s-s1)]:

zip(int,M,[s,s0,s1]);

Or, if you really do have Matrices (instead of lists as above),

M:=Matrix([[(s-s1)*(s-s2), (s-s0)*(s-s2), (s-s0)*(s-s1)]]):

zip(int,M,Matrix([[s,s0,s1]]));

acer

@Christopher2222 It might be more sensible to first ask whether it can be done, and save until later the asking as to how.

No, you cannot really add another subpackage (like Astronomy) "under" ScientificConstants. Well, not without unprotecting it and replacing & resaving the whole package with a modified version. And you wrote that you didn't want to do that. And it's not easiest, anyway.

It is not true that a separate package cannot make use of ScientificConstants routines, which you supposed in a comment responding to another answer.

The AddConstant help-page suggests placing definitions for new constants inside an initialization file. But that doesn't really bundle them nicely, for re-use by other people. You could, as an alternative, put such definitions (via AddConstant) into the ModuleLoad of a package which you savelib to a Library archive. Once you do that then whenever that package is first read by Maple in a session the definitions will be loaded.

The module below is not part (or a subpackage) of ScientificConstants.

First let's write the source of the module.

Astronomy:=module() option package;
export GetConstant,GetValue,GetError,GetUnit,Constant;
local ModuleLoad;
   ModuleLoad:=proc()
      ScientificConstants:-AddConstant(':-mass_of_Jupiter',
                                       'symbol'=':-MJupiter',
                                       'value'=1.90e27,'units'=':-kg')
   end proc:
   Constant:=ScientificConstants:-Constant;
   GetConstant:=ScientificConstants:-GetConstant;
   GetValue:=ScientificConstants:-GetValue;
   GetError:=ScientificConstants:-GetError;
   GetUnit:=ScientificConstants:-GetUnit;
end module:

Now let's savelib it to a Library archive.

try
  LibraryTools:-Create("c:/temp/astro.mla",WRITE):
catch "there is already an archive":
end try;

libname:="c:/temp/astro.mla",libname:
savelib(Astronomy);

Now let's restart, and load it from that archive.

restart:

libname:="c:/temp/astro.mla",libname:

with(Astronomy);

      [Constant, GetConstant, GetError, GetUnit, GetValue]

GetConstant(MJupiter);
                                                       27  
    mass_of_Jupiter, symbol = MJupiter, value = 1.90 10  , 

      uncertainty = undefined, units = kg

GetValue(Constant(MJupiter));

                                  27
                           1.90 10  

Note that the exports of the module above are not strictly needed. You could just as well load ScientificConstants and use its exports instead. I put in those exports just to save a step. The key bit is the ModuleLoad.

Notice that `evalf/Sum` fails here, and plot is going to try using `evalhf` and `evalf`, but not `value`.

I don't see what the problem is with `add`, except of course that it is not fastest. And you have to be careful about avoiding evalhf mode.

restart;

S := x->1+Sum(floor((x/(Sum(floor(cos(Pi*(factorial(n-1)+1)/n)^2),
              n = 1 .. r)))^(1/x)), r = 1 .. 2^x):

seq(value(S(x)), x = 1 .. 7);

                     2, 3, 5, 7, 11, 13, 17

seq(evalf(S(x)), x = 1 .. 7);

                  2., 3., 3., 5., 11., 13., 6.

S := x->1+add(floor((x/(add(floor(cos(Pi*(factorial(n-1)+1)/n)^2),
              n = 1 .. r)))^(1/x)), r = 1 .. 2^x):

seq(S(x), x = 1 .. 7);

                     2, 3, 5, 7, 11, 13, 17

As to the speed of this, you could try and optimize in at least two ways. The inner addition (in n) might be replaced with calls to a procedure that had `option remember`. And the outer addition (in r) could be done using a recursive procedure. This could help you avoid reproducing all previous work in order to generate each successive result.

Also, you should try and plot `S` not `S(x)`, since S is not prepared for symbolic `x` as argument (when using `add`). And, did you want to do a point-plot? I mean, do you really want to have this try and plot for non-posint x? Note that `evalhf` mode of floating point evaluation goes wrong for the `add` example I gave. But that is easily avoided.

plot(evalf@S,sample=[seq(i,i=1..7)],adaptive=false,
     view=[1..7,1..20]);

acer

Are you trying to use commands from the Physics package (eg. KroneckerDelta) but might have forgtotten to load that package?

Try issuing the command,

  with(Physics):

before your other work.

acer

How accurate do you need the approximation to be? Maybe you could consider numeric approximation instead of series.

Does it need to have a relative error of 1e-10 say? Or is it just for plotting, in which case an error like 1e-3 may be adequate? For the higher accuracy you could try the numapprox package and one of its routines such as minimax or chebyshev. For just plotting you could compute an Array of values and then interpolate with the CurveFitting package or other means.

Just how expensive  is to to do a single evaluation at a given value of x?

Could you not upload a worksheet containing the function, and attach it to your post?

acer

If you want to create an image format file of the plot, then start with the help-page for the plotsetup command.

If you want to export the data of a plot then you could look at plots:-getdata (or directly compute an Array of data yourself) and an exporting command such as ExportMatrix. For example,

f:=sin(x)+cos(Pi*x):

P:=plot(f,adaptive=false,sample=[seq(x,x=0..10,10/(5000-1))]):

plottools:-getdata(P)[3];
                          [ 5000 x 2 Matrix      ]
                          [ Data Type: float[8]  ]
                          [ Storage: rectangular ]
                          [ Order: Fortran_order ]

...and then you can apply the ExportMatrix (or writedata, etc) to that 5000x2 Matrix.

acer

If you know that the function is periodic, and if you know the period, then you should only need to evaluate the function enough to compute plot data over a single period. A compound plot over several periods can then be done without any further evaluation of the function, merely by a translation of the single-period plot.

pforce:=p*sin(omega1*t)+3*cos(3*omega1*t);

d:=[p=3,omega1=4.2];

period:=eval(2*Pi/omega1,d);

P:=plot( eval(pforce,d), t=eval(0..2*Pi/omega1,d) ):

P;

plots:-display( seq(plottools:-translate(P,period*i,0), i=-3..3) );

acer

There is a distinction between trying to cast the whole Matrix when it is needed as datatype=sfloat (simply use evalf on the no-datatype Matrix, and create a new Matrix), and casting each entry as it written into the Matrix. The former is likely more efficient if you only have to do it once but is costly if you have to keep repeating the action, while the latter adds a cost to every individual entry access. Your Question asks for the latter, it seems.

You could use a customized (user-defined) rtable indexing function as a corrective gatekeeper.

> restart:

> `index/coerce_float` := proc(idx,M,val)
> local g,idx1;
> idx1:=op(idx);
>     if nargs = 2 then
>       # retrieval
>       M[idx1];
>     else
>       g := op(val);
>       # storage
>       if type(g,sfloat) then
>           M[idx1] := g;
>       elif type(g,{hfloat,integer}) then
>           M[idx1] := SFloat(g);
>       elif type(g,realcons) then
>        M[idx1] := evalf(g);
>       else
>          M[idx1] := g; # watch it burn
>     end if;
>   end if;
> end proc:

> M := Matrix(2,2,'shape'='coerce_float','datatype'='sfloat');

                                     [0.  0.]
                                M := [      ]
                                     [0.  0.]

> M[1,1]:=-4.5:
> M;
                                 [-4.5  0.]
                                 [        ]
                                 [  0.  0.]

> lprint(M);
Matrix(2, 2, {(1, 1) = -4.5}, datatype = sfloat, storage = rectangular,
order = Fortran_order, shape = [coerce_float])

> M[1,1]:=sin(3):
> M;

                             [0.1411200081  0.]
                             [                ]
                             [          0.  0.]

> lprint(M);
Matrix(2, 2, {(1, 1) = .1411200081}, datatype = sfloat, storage = rectangular,
order = Fortran_order, shape = [coerce_float])

> M[1,1]:=HFloat(2.3):
> M;

                          [2.29999999999999982  0.]
                          [                       ]
                          [                 0.  0.]

> lprint(M);
Matrix(2, 2, {(1, 1) = 2.29999999999999982}, datatype = sfloat, storage = rectangular,
order = Fortran_order, shape = [coerce_float])

> M[1,1]:=Int(x^2,x=0..1):
> M;

                             [0.3333333333  0.]
                             [                ]
                             [          0.  0.]

> lprint(M);
Matrix(2, 2, {(1, 1) = .3333333333}, datatype = sfloat, storage = rectangular,
order = Fortran_order, shape = [coerce_float])

> M[1,1]:=x:
Error, (in index/coerce_float) unable to store 'x' when datatype=sfloat

acer

One easy way is to use "multiple assignment".  For example,

  evals, evecs := Eigenvectors(DD);

Another way is to use indexing.

  result := Eigenvectors(DD);

  result[2];

You should not be using `D` as a variable, especially at the top level. It is the name of the differential operator (Ie. it differentiates procedures.) That is why it is protected at the top level.

acer