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.

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

Unlike some of the other dials, the code inside DialKp1 does not have the line,

   kp1 := Do(%DialKp1):

When I add that line, the coefficient in the formula shown in MathBase gets updated as the dial is turned.

You should probably think about cleaning up the code in the Components. One easy improvement would be to get rid of the `restart`s and `with`s. You could also put code repeated amongst different sliders into centralized location(s) such as a procedure, a simple call to which would then be the only code inside each Component. You could put such procedure definitions inside a Code Edit Region, or the Startup code of the Document. Things might be easier for you, if you only had to make edits to one location instead of to all those different locations.

acer

Is this a kind of repeat of your previous question? Ie. the field directions for C vs I may not be constant and meaningful unless (possibly all) other variables in the system are constant. But since the solution to a differential equation system doesn't in general have all-but-two variables constant then there's no useful meaning to such a 2D planar slice of the field. (Anyone agree with me, or am I spouting nonsense?)

Note that Maple uses `I` to denote the square root of -1. You should choose another variable name, instead of trying to use that as a variable. (It is possible to adjust what maple uses to represent the imaginary unit, but as a new user of Maple I suggest instead keeping things simple and merely choose another name for your variable.)

acer

You're welcome for the answer. But you don't have to branch off the topic, just to add a followup comment. (Use the Reply, not the Branch link, below a post.)

I'm sure the biology behnd it is interesting. But it doesn't matter what the axes are called: x,y, and z, or healthy, infected, and immune. Those are just labels.

The field direction for healthy cells versus infected cells changes with the number of immune cells. There is no unique representation (in a plane representing healthy cells and infected cells only) for the field, given your original system of differential equations in all three quantities. The simple 2 dimensional answer that you're asking for, does not exist as far as I can see. Please see my earlier response.

acer

It's an interesting question: how to project the 3D spacecurve onto the x-y plane along with the field arrows corresponding to 3D points (only) along that curve. I'm not sure how to do that easily, without just building and inserting a bunch of arrows by hand.

The 3D field arrows corresponding to 3D points along a given vertical line (x and y both constant) do not all point in a constant x-y direction. So there is no obvious correct choice to make for the field arrows elsewhere in the x-y plane (except maybe for the x-y points of the projection of the 3D spacecure).

Here's the 3D plot. Perhaps it will convey the message you hoped for from the 2D portrait.

restart:

DE:={D(x)(t)=y(t)-z(t),D(y)(t)=z(t)-x(t),D(z)(t)=x(t)-y(t)}:

dsn:=dsolve(DE union {x(0)=1,y(0)=0,z(0)=2}, numeric, output=listprocedure):

plots:-display(
  plots:-odeplot(dsn, [x(t),y(t),z(t)], t=-2..2),
  plots:-fieldplot3d(subs(x(t)=x,y(t)=y,z(t)=z,
                        eval([diff(x(t),t),diff(y(t),t),diff(z(t),t)],
                             convert(DE,diff))),
                   x=-0.5..2.5,y=-0.5..2.5,z=-1.0..3.0) ,axes=box);

acer

As given, the worksheet attempts to do symbolic solving of the DE. This results, after using quite a lot of time and memory, in a very long solution which could be of less practical use than you might want. (Trying to do symbolic solving of equations with floating-point coefficients is often a warning that something is amiss, IMO.)

You could try something like this, at the end,

sol3 := dsolve([eqB, eqC, eqD, eqE, eqF,
                speed(0) = 0, V(0) = 0.7047824398e-1],
               [speed(t), Tm(t), Tw(t), V(t), omega(t)],
               numeric, output = listprocedure):

spd := subs(sol3, speed(t)):

pos := proc(T) evalf(Int(spd,0..T,epsilon=1e-3)); end proc:

st:=time():
plot([spd,pos], 0..20);
time()-st;

That's computing both the speed and position numerically. Now, didn't Markiyan recently answer a related earlier question by you, citing an old post in which Joe Riel suggested augmenting the DE system so as to allow dsolve/numeric itself to do the integration (as opposed to using evalf/Int to do it with numeric quadrature)? That was probably good advice, because dsolve/numeric is supposed to be clever about generating a sequence of computed values, whereas the above scheme is going to repeat effort when integrating numerically from T=0 out to T=blah, then once more from T=0 to further point T=blech, etc, etc for all the points needed for the plot.

If you choose to stick with a numeric scheme for calculating position like I show here, then Patrick T's related advice (also to you) seems to save about 1/4 the computational effort while plotting (but this may involve luck, that the repeated numeric quadratures re-use points in time -- instead of using all ugly unre-usable gaussian points, say...). For example, to compare, this seems to take 3/4 the time of the above plot,

spd := subsop(3=remember, subs(sol3, speed(t))):
st:=time():
plot([spd,pos], 0..20);
time()-st;

Your attached worksheet seems to contain some mysterious hidden call which starts the ODE assistant (maplet), if run using the triple-exclam from the menubar. Weird. Here's an edited version using what I've mentioned. (I didn't hunt for the mysterious thing.)

Do you still need to resolve the data table interpolation question, afer trying this? On-the-fly interpolation shouldn't be faster than your degree 6 polynomial fit. The only practical benefit I can think of for still trying it is if you feel it would interpolate the data better. Less wiggle, or what have you.

Gear_Ratio_Study-ne.mw

acer

Yes.

You can use procedures for this. Put all the computational code inside a procedure whose parameters are the variables you wish to vary. Then, instead of assigning some values to those variables at the top-level (ie. instead of outside the proc), call your procedure with those values. Or call it again with different values. And so on.

But it's not a trick! It's the natural progression from the basic assign-values-then-loop-while-doing-stuff-then-plot-and-rely-on-menubar-triple-exclam-to-execute-all approach to using Maple.

Sure, you might have to cut some of the material and paste it into just one execution block to contain the procedure.

If you run into problems in accomplishing this, you could upload the worksheet in a new comment to this thread and tell us what you want to vary.

acer

Try command-completion, while typing. That also works for me in Maple 12 (which I believe you still use.)

In other words, do it like so:

• type in diff
• hit the command-completion key sequence (Ctl-space on Windows)
• Choose either the third or second entry in the list the pops up, ie. diff (inline), or diff (inline partial)
• Don't put the expression to be differentiated in the numerator of the "derivative symbol", put it afterward.

You want to end up with some like

and not something like

If you want to control how much of an expression gets the applied to it, then use round brackets to delimit the scope. Note the parsed meaning,

rather than being equal to,

You can also use the appropriate entry from the Expression palette.

acer

DEtools[DEplot] is really just out of date and inefficient in its technique.

A leaner and faster approach is to use plots:-fieldplot alongside a sequence of (with varying ICs) of plots:-odeplot calls.

I'm going to branch off this Answer, and show both plots:-animate and Embedded Component implementations of both the old DEplot and that new "combined" approach.

As for saving the worksheet and huge memory use by the GUI, well you can solve that by removing all the output from the worksheet before ever saving it. You might object and say, "but my animation takes 45 sec to create on the very fastest machine! I don't want to have to wait every time I open it!" To which I'd respond, just look how much faster the combined approach is: maybe 40-60 times faster to create the whole animation, for the ODE system in this thread.

And the Embedded Component implementation saves, with output included, with no problem. And the combined approach allows it to display smoothly and quickly, while allowing for nice manipulation of two parameters.

acer