The data that you generate (x and f(x)) will have to be in a data structure in order to get plotted. You can insert them into a preconstructed Matrix, or form them as a list of lists, or even print and subsequently read them (into a Matrix) from a file.

Here are a few ways, (trying to keep to your loop methodology where appropriate and practical),

restart:

M:=Matrix(10001,2):
for x from 0 to 100 by 0.01 do
   M[trunc(x*100)+1,..] := <x,sin(x)>;
od:
plot(M);

N:=Matrix(10001,2):
for x from 0 to 100 by 0.01 do
   idx := trunc(x*100)+1;
   N[idx,1], N[idx,2] := x, sin(x);
od:
plot(N);

P:=Matrix(10001,2):
for i from 1 to 10001 do
   x := (i-1)*0.01;
   P[i,1], P[i,2] := x, sin(x);
od:
plot(P);

L:=[seq([x,sin(x)],x=0..100,0.01)]:
plot(L);

#FileTools:-Remove("thedatafile");
for x from 0 to 100 by 0.01 do
  fprintf("thedatafile", "%.10f, %.10f\n",x, sin(x));
od:
fclose("thedatafile");
Q := ImportMatrix("thedatafile", source=csv):
plot(Q);

LinearAlgebra:-Norm(M-Matrix(L));
LinearAlgebra:-Norm(M-N);
LinearAlgebra:-Norm(M-P);
LinearAlgebra:-Norm(M-Q);

acer

Since the upper limit of your integral will evaluate to a float (for nice numeric delta) then `int` will call out to do numeric quadrature. But better to keep things straight and use `Int`, so as to ensure that `int` never wastes time trying to do a mixed symbolic-float integral in some "exact" way (which is often a Bad Idea, as well as time wasting when it fails).

Also, NextZero may try an bump up Digits, which can be unfortunate for evalf(Int(...)) if it makes Digits too high for fast double precision computation, or it this interferes with the convergence acceptance. So I use a proc `fdf` which forces Digits=15 each time its called.

I find that these times below can vary from run to run, even when trying to measure the timings while clearing remember tables. (I suspect there may be a remember table that I am not always clearing.)

f:=Int((50-50*delta)*(50+50*delta)*(0.3e-1+(.36*(1-exp(-.39*t)))
       /(12.00000000*exp(-.39*t)+1))
       *(1-(1-exp(-.39*t))/(12.00000000*exp(-.39*t)+1))/1.08^t,
       t = 0 .. -2.564102564*ln((1.+5.*delta)/(53.+5.*delta)))
     +(500.*(-1.-5.*delta))*(-(1.*(1.+4.*ga))/(-53.+48.*ga))^.197336002905
      /((-1.+1.*ga)*(1.+delta)*((1.+5.*delta)/(53.+5.*delta))^.197336002905):

#plots:-display(seq(plot(eval(f,ga=m),delta=0..2),m=0.1..0.8,0.1));

df:=diff(eval(f,ga=0.1),delta):

forget(evalf); forget(`evalf/int`):st:=time():
fsolve(df=0,delta=0..1);
                          0.1706030333
time()-st;
                             0.219

fdf:=proc(m) Digits:=15: eval(df,delta=m); end proc:

forget(evalf); forget(`evalf/int`):st:=time():
RootFinding:-NextZero(fdf,0.0);
                          0.1706030333
time()-st;
                             0.156

forget(evalf); forget(`evalf/int`):st:=time():
fsolve(fdf,0..1);
                          0.1706030333
time()-st;
                             0.187

A little more time (maybe 15% off those times above) could be saved by specifying some evalf/Int options (such as method=_d01ajc, and maybe digits=15 and epsilon=1e-13). But it's awkward to get these optiona into the Int() call, since they muck with the `diff` call if they are present in the original `f`.

acer

I don't know whether this will help in your more complicated examples...

> restart:

> S:=sum(sum((b[i]-b[j]),i=1..n),j=1..n);

                            n   /          /  n       \\
                          ----- |          |-----     ||
                           \    |          | \        ||
                            )   |          |  )       ||
                     S :=  /    |-n b[j] + | /    b[i]||
                          ----- |          |-----     ||
                          j = 1 \          \i = 1     //

> with(student):

> expand(subs(sum=Sum,S));

                     /  n       \   /  n       \ /  n    \
                     |-----     |   |-----     | |-----  |
                     | \        |   | \        | | \     |
                     |  )       |   |  )       | |  )    |
                  -n | /    b[j]| + | /    b[i]| | /    1|
                     |-----     |   |-----     | |-----  |
                     \j = 1     /   \i = 1     / \j = 1  /

> simplify(value(expand(subs(sum=Sum,S))));

                                      0

acer

The single-backslash \ is being interpreted as an escape character. Try double-backslash like \\, or forward slash.

acer

I suspect that what is happening is that some (or at least one) of the calls to fsolve is not succeeding. When it fails to find a root, it returns as an unevaluated call (to itself). This is not numeric.

It only takes one non-numeric value in the Matrix for `matrixplot` to complain. (It doesn't matter whether some values of the Matrix, such as the [1,1] entry, are numeric. The `matrixplot` command expects all entries to be numeric. So even just one failing fsolve attempt in your code -- as it's currently written -- is going to ruin the `matrixplot` attempt.

You could put a test into your code, to check that Temp is of the expected type (eg, type(Temp,identical(T)=numeric) or what have you). See here (and esp. Joe's comment) for a note on that kind of check outside of a procedure. If any fsolve calls fail then you could try another root-finding mechanism, or insert a special value into the Matrix, or print a Warning, or write nothing new to the relevent Matrix entry, etc. You could code that kind of fallback action with a conditional on the type-check.

Could you please try not to post multiple Questions for the same essential issue? Thanks. If you have followup details or comments on a single issue (at essence) then you can just add a Comment to the original.

acer

The long ton is already present.

> restart:

> kernelopts(version);

          Maple 12.02, X86 64 WINDOWS, Dec 10 2008 Build ID 377066

> convert(1,units,ton[long],lb);

                                    2240
> restart:

> with(Units:-Standard):
> Units:-UseSystem(FPS);

> 2200*Unit(lb/s)+40*Unit(lb/s);

                                         /lb\
                                2240 Unit|--|
                                         \s /

> Units:-AddSystem(mysys,Units:-GetSystem(FPS),ton[long]);
> Units:-UseSystem(mysys);

> 2200*Unit(lb/s)+40*Unit(lb/s);

                                   /ton[long]\
                               Unit|---------|
                                   \ second  /

acer

I've been investigating this idea in some spare time, for a little while. I believe that I, be found a way. At the same time, I've been looking at animatio in Plot Components (with multiplie sliders) or using the "old" animate mechanisms. It's looking promising, although there are some GUI memory leaks that need fixing (but the need fing for plain old animat of 3d plot animate too!). But I'm on vacation for the next 8 days, without a computer or access. Here is a hint, for any adventurous code who has studied theExplore Library code: instead of loading the XML of the new exploring components into a new worksheet, load into a hdb to replace a pre-existing Task,xxxx topic, and then open that topic. Almost all the steps except Insert default Content can be programmed to happen without ned for mouse interaction!d

acer

Have you looked at (right-click) context-sentive menus, with annotations enabled as a GUI option?

One can customize those menus. It might be possible to kludge a set of new/modified context-menu actions which have close to what you want in look and feel. (Possibly with some print-extension stuff.)

Describing to us a complete and explicit example, from the simple algebra exercises say, would help.

Here's the kind of issue that might arise, and it'd be key to know how you'd like it handled: Will you need to simultanelously add and subtract the very same term to the RHS, where the added term is combined with the RHS and the subtracted term is left to augment the new resulting RHS? Do you need both to appear simultanteously on one line, before any combination? If so, then special means have to be used so that they can both be printed ( the `x` and the `-x`) without immediately cancelling. Eg,  LHS = RHS + x - x   as displayed output.

If I understand, you don't want any Maple commands around. Just high school math. So the "steps" returned by the Equation Manipulator assistant is likely not going to help you.

acer

It looks like the more useful thing for you to do would be to remove that assignment to wm.

 (**) restart:

 (**) eq:=diff(theta(t),t)*y(theta(t))+x(t)=sin(x(t));

                    /d          \
              eq := |-- theta(t)| y(theta(t)) + x(t) = sin(x(t))
                    \dt         /

 (**) wmeq:=wm=diff(theta(t),t);                               

                                        d
                           wmeq := wm = -- theta(t)
                                        dt

 (**) subs((rhs=lhs)(wmeq),eq);    

                       wm y(theta(t)) + x(t) = sin(x(t))

 (**) subs(wmeq, %);
                 /d          \
                 |-- theta(t)| y(theta(t)) + x(t) = sin(x(t))
                 \dt         /

But maybe aliasing is something else you could consider.

 (**) restart:

 (**) eq:=diff(theta(t),t)*y(theta(t))+x(t)=sin(x(t));

                    /d          \
              eq := |-- theta(t)| y(theta(t)) + x(t) = sin(x(t))
                    \dt         /

 (**) alias(wm=diff(theta(t),t)):

 (**) eq;

                       wm y(theta(t)) + x(t) = sin(x(t))

 (**) diff(eq, t);

/ 2          \
|d           |                 2                  /d      \
|--- theta(t)| y(theta(t)) + wm  D(y)(theta(t)) + |-- x(t)| =
|  2         |                                    \dt     /
\dt          /

              /d      \
    cos(x(t)) |-- x(t)|
              \dt     /

It depends on what end you're trying to attain.

acer

kernelopts(printbytes=false)

It's slightly harder to find, because for some reason it's a kernelopts thing and not an interface thing.

acer

No, Maple isn't linked to the ACML.

 

(At some point, the Maple-Nag Connector toolbox on 32bit Linux allowed the choice of using a version of the NAG C Library Mark 8 linked against the ACML. But Maple itself wasn't/isn't linked against it, not even from the portions of NAG which have been more seamlessly incorporated for use by Library and kernel.)

acer

You can do this with a PlotComponent.

Here is a simple example. There's a lot of choice in how to do it. (ie. Toggle buttons, buttons which change their own captions between "Play" and "Pause" and do multiple duty, etc, etc.)

Right-click on the button and the checkbox to see their underlying code (Action when Clicked) which in both cases are just simple calls to the two routines RunMe and PauseMe.

 

>  restart: >  # You only need to set up these procedures once. # You can hide this in your StartUp region, or in a # Code Edit Region. RunMe:=proc(plotcomp, checkcomp) uses DocumentTools;    SetProperty(checkcomp,'value',"false");    SetProperty(plotcomp,'play',"true"); end proc: PauseMe:=proc(plotcomp, checkcomp) uses DocumentTools;    if GetProperty(checkcomp,'value')="true" then       SetProperty(plotcomp,'play',"false",'refresh'=true);    else       SetProperty(plotcomp,'play',"true",'refresh'=true);    end if; end proc: >  # Create the animation. anim:=plots:-animate(plot,[sin(x*a),x=-20..20], a=-10..10): >  # Insert the animation into a Plot Component. DocumentTools:-SetProperty('Plot0','value',anim,'refresh'=true); >  >    >  RunMe('Plot0','CheckBox0'); # runs it programmatically >

 

Download component_anim.mwcomponent_anim.mw

acer

Just because `evalf` itself is a builtin doesn't mean that it cannot call back to interpreted Library level routines.

showstat(`evalf/Sum`);

showstat(`evalf/Sum1`);

showstat(`evalf/Sum/LevinU_Symb/1_inf`);

and so on..

acer

You say you would like the plot "accurate", but could you quantify this?

You can adjust the tolerance (option `epsilon`) below. I got about a factor of 2 speed up by simplifying the integrand. Another 3-fold speed up or so came from relaxing the quadrature tolerance (it still seems pretty accurate) and by specifying the method. Then, since a sample of 6 data points seem to show a nice, moderately gentle, monotonic function I generated an interpolated smooth plot, which took less than 10 seconds on a fast machine.

The bottom line cost (not counting the 4 sample points for speed comparison) was roughly: 8 to 9 seconds to get the plot, after about a second to generate the equations do the symbolics.

restart:

myH1:=(x,n,t)->HankelH1((1/2)*n-1/2,x*exp(-t)):
myH2:=(x,n,t)->HankelH2((1/2)*n-1/2,x*exp(-t)):
myHP1:=(x,n,t)->diff(myH1(x,n,t),x):
myHP2:=(x, n, t)->diff(myH2(x, n, t), x):
myBeta:=(x, n, t)->((1/2)*n-1/2-I*x)*myH1(x, n, t)+x*myHP1(x, n, t):
myBetastar:=(x, n, t)->((1/2)*n-1/2+I*x)*myH2(x, n, t)+x*myHP2(x, n, t):
BE:=(x, n)->x^(n-1)/(exp(x)-1):
ParticleDensity:=(x, n, t)->BE(x, n)*myBeta(x, n, t)*myBetastar(x, n, t):

ParticleNumber:=t->Re(evalf(Int(ParticleDensity(x,3,t),x=0..infinity))):

forget(evalf): forget(`evalf/int`):
st:=time():
ParticleNumber(0);

                          0.1863252512

time()-st;

                             8.612

ParticleNumberExpr := ParticleDensity(x, 3, t):
ParticleNumberExpr := simplify(ParticleNumberExpr) assuming t>0, t<1, x>0, x<infinity:
ParticleNumber2:=subs(GG=Int(ParticleNumberExpr,x=0..infinity,
                             epsilon=1e-5,method=_CCquad),t->Re(evalf(GG))):

forget(evalf): forget(`evalf/int`):
st:=time():
ParticleNumber2(0);

                          0.1863252492

time()-st;

                             1.092

forget(evalf): forget(`evalf/int`):
st:=time():
ParticleNumber2(1); # now the alternate approach first

                          7.773916361

time()-st;

                             1.450

forget(evalf): forget(`evalf/int`):
st:=time():
ParticleNumber(1);

                          7.773916460

time()-st;

                             9.204

st:=time():
Tin:=Vector(6,(i)->(i-1)/5,datatype=float[8]):
Tout:=Vector(6,(i)->ParticleNumber2((i-1)/5),datatype=float[8]):
#plot(); # without interpolation, just a point-plot
Func:=x->CurveFitting:-ArrayInterpolation(Tin,Tout,x,method=spline):
plot(Func,0..1);

time()-st;

                             8.627

I haven't looked to whether you should/could be computing Re of the numeric integral, or the numeric integral of Re of something, etc. And the heavy use of operators, for everything that could just as well be an expression, is not to my taste. But I figure you know what you want, and that that's ok.

acer

The Roots command uses fsolve and its `avoid` option internally to do its work.

The methods fsolve will bring to bear are mostly all iterative schemes whose success depend on the chosen initial point. And, internally, fsolve may use some finite number of initial starting points.

But when the specified candidate range gets large, and when there are roots close enough together, then it becomes possible that for a particular root none of the generated initial starting points will converge to it. This can be true for a set (uniform, or other) collection of initial points, and it can also become likely if the initial points are generated randomly.

In other words, for a large interval in the domain the relative size of a basin of convergence for some particular root might be relatively very small. And with a fixed number of generated initial points it may be that (few or) none converge to the given root.

Hence, one possible rule to try and follow would be to keep the stated range (in the variable, in the domain) as short as reasonably possible.

An alternative is to use RootFinding:-NextZero, which searches for the very next root while moving to the right (in the direction +infinity on the real axis).

I think that the NextZero routine has improved in recent releases. Maybe it should be considered for use in Student:-Calculus1:-Roots.

Here is a simple procedure which uses NextZero to try and find all real roots of expression `expr` within the real range `a` to `b`. When it works, it seems to be quite a bit faster than Roots.

restart:

findroots:=proc(expr,a,b,{guard::posint:=5,maxtries::posint:=50})
local F,x,sols,i,res,start,t;
   x:=indets(expr,name) minus {constants};
   if nops(x)>1 then error "too many indeterminates"; end if;
   F:=subs(__F=unapply(expr,x[1]),__G=guard,proc(t)
      Digits:=Digits+__G;
      __F(t);
   end proc);
   sols,i,start:=table([]),0,a;
   to maxtries do
      i:=i+1;
      res:=RootFinding:-NextZero(F,start,
                                 'maxdistance'=b-start);
      if type(res,numeric) then
         sols[i]:=fnormal(res);
         if sols[i]=sols[i-1] then
            start:=sols[i]+1.0*10^(-Digits);
            i:=i-1;
         else
            start:=sols[i];
         end if;
      else
         break;
      end if;
   end do;
   op({entries(sols,'nolist')});
end proc:

findroots(x^3-3*2^x+3, -1000, 1000);

          -1.189979747, -0., 2.226983010, 6.592619383

[edited: I edited the maxdistance value in the call to NextZero above, to make it simpler and to better handle the first time in the loop.]

Almost any practical numerical rootfinding scheme can be broken. Hence, while suggestions for improvements are welcome, there's no bounty for examples which break this scheme. With that in mind, here is a performance example. The first columns is the total number of roots, the second column is the number of roots found, the third column is the maximal forward error by resubstituting the roots into the expression, and the fourth column is the elapsed time.

V:=LinearAlgebra:-RandomVector[row](100,generator=-1.0..1.0):

for i from 1 to 10 do
  expr:=expand( mul((4.3)^x-(4.3)^V[j],j=1..i) );
  d1:='d1':d2:='d2':d3:='d3':gc():
  st:=time();
  sol:=Student:-Calculus1:-Roots(expr,x=-1..1,numeric);
  maxerr:=max(seq(abs(eval(expr,x=X)),X=sol));
  print([i,nops(sol),nprintf("%.1e",evalf[2](maxerr)),time()-st]);
end do:

                     [1, 1, 0.0e+00, 0.094]
                     [2, 2, 1.0e-10, 0.203]
                     [3, 3, 2.0e-10, 1.107]
                     [4, 4, 4.0e-10, 1.638]
                     [5, 5, 7.0e-10, 1.186]
                     [6, 6, 1.4e-08, 1.872]
                     [7, 7, 8.7e-06, 2.449]
                     [8, 8, 5.2e-06, 4.352]
                     [9, 9, 1.5e-04, 5.008]
                    [10, 10, 5.9e-05, 7.425]

for i from 1 to 10 do
  expr:=expand( mul((4.3)^x-(4.3)^V[j],j=1..i) );
  d1:='d1':d2:='d2':d3:='d3':gc():
  st:=time();
  sol:=findroots(expr,-1,1,guard=7);
  maxerr:=max(seq(abs(eval(expr,x=X)),X=[sol])):
  print([i,nops([sol]),nprintf("%.1e",evalf[2](maxerr)),time()-st]);
end do:

                      [1, 1, 1.0e-10, 0.]
                      [2, 2, 0.0e+00, 0.]
                      [3, 3, 2.0e-10, 0.]
                      [4, 4, 3.0e-10, 0.]
                     [5, 5, 7.0e-10, 0.016]
                     [6, 6, 4.1e-07, 0.015]
                      [7, 7, 5.1e-07, 0.]
                      [8, 8, 7.5e-06, 0.]
                      [9, 9, 1.4e-04, 0.]
                    [10, 10, 5.9e-05, 0.016]

NB. It's quite possible that the findroots procedure is not as careful as it might be about avoiding a runaway situation wherein its main loop stops only due to the maxtries parameter.

acer