Did you not already ask that here, and see the posted answer?

nb. As mentioned there, a list in Maple is entered using [] square brackets, not {} brackets which are for sets.

acer

The HFloat underperforms compared, say, to immediate integers.

restart:
kernelopts(gcfreq=[10^7,0.0000001]):
st:=time():
for i from 1 to 5 do
  for j from 1 to 10^6 do
    HFloat(j);
  end do;
  gc();
  print([kernelopts(bytesalloc),kernelopts(bytesused),kernelopts(gctimes)]);
end do:
time()-st;
                  [104903924, 4037634432, 341]
                  [104903924, 4076101300, 342]
                  [104903924, 4113039968, 343]
                  [104903924, 4150709568, 344]
                  [104903924, 4188603216, 345]
                             7.515

restart:
kernelopts(gcfreq=[10^7,0.0000001]):
st:=time():
for i from 1 to 5 do
  for j from 1 to 10^6 do
    j;
  end do;
  gc();
  print([kernelopts(bytesalloc),kernelopts(bytesused),kernelopts(gctimes)]);
end do:
time()-st;
                   [851812, 4192460624, 347]
                   [851812, 4192462656, 348]
                   [851812, 4192464632, 349]
                   [851812, 4192466544, 350]
                   [851812, 4192468488, 351]
                             0.719

restart:
kernelopts(gcfreq=[10^7,0.0000001]):
st:=time():
for i from 1 to 5 do
  for j from kernelopts(maximmediate) to kernelopts(maximmediate)+10^6 do
    j;
  end do;
  gc();
  print([kernelopts(bytesalloc),kernelopts(bytesused),kernelopts(gctimes)]);
end do:
time()-st;
                  [61723608, 4257649624, 353]
                  [61723608, 4285651840, 354]
                  [61723608, 4313653832, 355]
                  [61723608, 4341655848, 356]
                  [61723608, 4369657836, 357]
                             3.422

acer

When you use the "Apply a Command" entry from the context-sensitive menus, you can use this command in the pop up entry box (typed just as you see it below, including brackets):

(evalc@abs,evalc@argument)

Did you know that you can also customize the context-menus, by adding (new, or more complicated) commands, optionally with their own "annotations above the arrow"? So if you want you can add a new item that does just the phase conversion you want, in one shot. And it could be annotated however you wanted. And you could put it in your Maple initialization file so that it was available in all you new sessions, or in a particular Document's "start up region" so that anyone else could re-execute it. See here for some notes on customizing the context menus. Ask for help, if you want to try it but run into difficulties.

acer

Perhaps your recursive binary search routine could pass 4 arguments instead of 3 arguments, and always pass along the full L. (I realize that this might not be the original assignment, but I'll lay out the idea anyway, for efficiency's sake.)

ZeFouilleDicho( L, x, G, D)

Here, the G is the left end position of the subrange of L on which one wishes to search. And D is the right end position.

One starts it off with an initial call like ZeFouilleDicho( L, 4, 1, nops(L)) since nops(L) is the last position of the previously sorted list L.

The routine determines which half of L, ie. the subrange L[G..D], one need to recurse into. Of course L[G..D] is never formed or passed along. Only the full list L is ever re-passed. That's the efficiency point. (One could also get rid of the recursion altogether, and do a binary search. But I'll retain the recursive nature, in the spirit of the assigned task.)

Note that you might consider using square brackets like L:=[1, 2, 4, 5, 6, 9, 11, 16] to make L a Maple list rather than squiggly brackets {} which are for creating Maple sets. That's because you might want to extend the code to handle exact symbolic constants like sqrt(2) or a mix of floats and exact integers. Compare,

> L:={1.2, 3, 4.6, 5, 7}; 
                      {3, 5, 7, 1.2, 4.6}
> L[1..3]; # oops, L isn't sorted as one might mistakenly expect
                           {3, 5, 7}

> L:=[1.2, 3, 4.6, 5, 7]; 
                      [1.2, 3, 4.6, 5, 7]
> L[1..3];
                         [1.2, 3, 4.6]

This recursive algorithm in which L is only ever passed in full is interesting because it avoids creating the sublists like L[G..D] explicitly. It passes along the extra 2 (G & D) pieces of positional information which denote that subrange. By avoiding creating each sublist explicitly, its avoid creating a whole bunch of sub-lists as collectible garbage. (This is the situation described in the 3rd paragraph here, where you yourself own and author `f`.)

Here's one possible implementation (of my suggestion, and technically not of what's asked above):

bs := proc(L,x,g,d)
  local mp;
  if d=g then d;
  elif x=L[g] then g;
  elif x=L[d] then d;
  elif x<L[g] or x>L[d] then
    error "x not in L";
  else
    mp := trunc((d+g)/2);
    if g=mp then g,d;
    elif x=L[mp] then mp;
    elif x>L[mp] then
      if d=mp+1 then mp,d;
      else
        procname(L,x,mp,d);
      end if;
    else
      if g=mp+1 then d,mp;
      else
        procname(L,x,g,mp);
      end if;
    end if;
  end if;
end proc:

L:=[3,4,5,6,7,9,10,11,12]:
bs(L,9,1,nops(L));
                               6
bs(L,5.5,1,nops(L));
                              3, 4
bs(L,3,1,nops(L));
                               1
bs(L,12,1,nops(L));
                               9

N := 10^7:
L := [seq(i,i=10^3..N)]:
st,ba,bu := time(),kernelopts(bytesalloc),kernelopts(bytesused):
bs(L,L[43*N/100],1,nops(L));
time()-st, kernelopts(bytesalloc)-ba ,kernelopts(bytesused)-bu;

                            4300000
                         0.172, 0, 7804

Note that bytesused goes up very little, indicating that little garbage collection is needed and done. (Which helps keep the timing down.) Of course a Compiled routine could be faster, but that's a whole other game. The performance of 'bs` is not too bad, compared to ListTools:-Search which uses the very fast kernel builtin `member` routine. (It likely could be made better, if the recursion mandate were dropped.)

st,ba,bu := time(),kernelopts(bytesalloc),kernelopts(bytesused):
ListTools:-Search(L[43*N/100],L);
time()-st, kernelopts(bytesalloc)-ba ,kernelopts(bytesused)-bu;
                            4300000
                         0.093, 0, 1880

acer

> M := sort([5, 3, 9, 48, 97, 7, 100, 88]):

> with(ListTools):

> L:=(M-Rotate(M,-1)):
> M[Search(max(L),L)-1];
                               48

Since some or all of the other suggestions create temporary lists (rather than just walking the sorted list), I figured that one more would be in the spirit of fun.

acer

While the problem has already been solved, and explained, I wanted to briefly mention an alternative. It merely allows for an alternate syntax in which there is no burden to think up and supply an intermetiate dummy name (y, or s, above), and in which a substitution formula (including its reversal) is only needed once.

> with(Student:-Precalculus):

> thaw(algsubs(x^2=freeze(x^2), 'CompleteSquare'(x^4+2*x^2, x^2) ));
                                 2    
                         / 2    \     
                         \x  + 1/  - 1

Naturally, some people will prefer one method, and some people another. And it can depend on the task at hand.

acer

Duncan's answer looks like it might be close. (djc's answer is also very nice, except that it may be too advanced and a little onerous for the "regular" or new user to go about constructing such `diff/f` extensions -- esp. a lot of them, programatically.)

Maybe this variant might also give the Poster some ideas. Do this below in 2D input mode, in a Document.

> f := t -> t^2:

> Diff(f(x), x) = value(Diff(f(x), x));

                          d  / 2\      
                         --- \x / = 2 x
                          dx           

When entering the above, after typing the keys for, "Diff", immediately hit Ctl-spacebar so as to get command-completion. (That's the keyboad acceleration on Windows. On Linux it is Ctl-shift-spacebar I think. Not sure about OSX.) The idea is to select the second menu choice for the command completion, which is the inline Diff. Doing it this way allows you to get the 2D Math input to appear in the prettyprinted typeset way (unlike what shows in this response of mine here). When the command-completion is accepted, the cursor is at the "x" field in d/dx. Type x in that field. Then use right-arrow to move the cursor out, and then enter the f(x). All quite natural and pretty easy, I think.

Of course, you can also do it as above, but without having assigned to `f`. To do it that way, it may work ok if done like (the 2D input and command-compeleted version of),

> Diff(f(x),x) = value(eval(Diff(f(x),x),f=(t->t^2)));

                          d            
                         --- f(x) = 2 x
                          dx           

And, naturally, you can use unapply(t^2,t) instead of (t->t^2), or make some big statement generator proc out of all that if you want it to accept t^2 as as programmatic argument. Etc, etc.

I don't think that (all) inert operators are easily or properly available in the palettes. Maybe Diff is, in some obscure palette. But command-completion is easier I think, because the approach is more uniform. You don't have to flit around looking for inert palette entries for Int, Diff, Sum, etc. (And they probably may not even be all hooked up properly, anyway.)

I very much prefer 2D Math command-completion over palette entries. And, for me, manual mouse-selection followed by context-menu conversion from 1D input to 2D input is least preferred as it entails both full typing and the most mouse work.

Several of the nicer (IMO, more useful) answers have involved using value(), and not actual assignment as was literally requested by the Poster. There is a way to do the actual assignent. That could be done by entering the 2D input of d/dx f(x), and then converting it to a Maple literal name using the context-menu action 2-D Math -> Convert To -> Atomic Identifier. I agree with the earlier responders, who seem to have implicity stamped this as the wrong approach (by virtue of suggesting alternatives). Also, one could get into some unattractive situations involving quotes, which just seems better avoided.

acer

For an explicit formula such as you have presented, in which there is a formula giving V in terms of the other variables, I would suggest that `eval` is a suitable tool for the task. Even it can be used in several slightly different ways to get the result. (I'll separate them into distinct sessions, separated by `restart`, to illustrate,

> restart:
> eqn := V = U + a*t;
                          V = U + a t

> eval(eqn, [U=3, a=9.81, t=2] );
                           V = 22.62

> restart;
> V := U + a*t;
                            U + a t

> eval(V, [U=3, a=9.81, t=2] );
                             22.62

For an explicit formula such as you've given, even a simple command like `eval` is not necessary, provided that you make assignments. Ie,

> restart:
> V := U + a*t;
                            U + a t
> U:=3:
> a:=9.81:
> t:=2:

> V;
                             22.62

Of course, for other formulations (with V mixed into the RHS say) there are Maple commands like `isolate` which can be used. Or you can hit it with `solve` which can be a powerful hammer. For formulations in which V is only given or attainable by an implicit formula, a purely numeric solver like `fsolve` can be a heavy hammer.

acer

Using `add` instead of `sum` makes a big difference here. For such adding of a finite (explicit) number of terms, you should be using `add`.

Suppressing printing of results by terminating statements with full colon : instead of semicolon ; can also help. (Not here, since the example below has length too long for display in the Standard GUI at default settings. But if the output is going to be too long for you to sensibly view, then don't print it out.)

> restart:
> f:=x->add(n*x*sin(n*x),n=1..1000):
> g:=taylor(f(x),x=1,60):
> kernelopts(bytesalloc);
                           287257216

That's about 300MB of memory allocated. It was reasonably fast, too.

acer

How about 219978219978219978219978219978

> x:=219978219978219978219978219978:
> length(x);
                               30

> 4*x - parse(StringTools:-Reverse(convert(x,string)));
                               0

Some patterns such as below can be used in a repetitive way.

> for N from 4 to 9 do
>   sol:=Optimization:-Minimize(
>     1,
>     {4*add(X[i]*10^(i-1),i=1..N)=add(X[N+1-i]*10^(i-1),i=1..N),
>     seq(X[i]<=9,i=1..N), X[1]>=1, X[N]>=1},
>     assume=nonnegint, depthlimit=1000, integertolerance=1e-15);
>   x:=add(rhs(sol[2][i])*10^(i-1),i=1..N);
>   print(x,  parse(StringTools:-Reverse(convert(x,string))));
> end do:

                           2178, 8712
                          21978, 87912
                         219978, 879912
                        2199978, 8799912
                       21782178, 87128712
                      217802178, 871208712

acer

This behaviour is documented by the option stoppingcriterion on the Bisection help-page. By default, a relative test between successive x-values is the (only) stopping criterion used. But a function-value test is an alternative stopping criterion that can be specified.

This is part of a Student package, so it is intended to allow the student some means to vary the behaviour and observe effects. That's why it provides alternate behaviour specified by such options.

The default behaviour is a relative test (only) using a default tolerance (both documented),

> restart:
> p:=expand( (x-3)^3 );
                      3      2            
                     x  - 9 x  + 27 x - 27

> Student:-NumericalAnalysis:-Bisection(
>              p, x=[0,12],
>              output=sequence); 

   [0., 12.], [0., 6.000000000], [0., 3.000000000], 

     [1.500000000, 3.000000000], [2.250000000, 3.000000000], 

     [2.625000000, 3.000000000], [2.812500000, 3.000000000], 

     [2.906250000, 3.000000000], [2.953125000, 3.000000000], 

     [2.976562500, 3.000000000], [2.988281250, 3.000000000]

And now, with a function-value stopping criterion,

> Student:-NumericalAnalysis:-Bisection(
>              p, x=[0,12],
>              output=sequence, stoppingcriterion=function_value); 

           [0., 12.], [0., 6.000000000], 3.000000000

Thus rests a defense of sorts. Now a counterpoint or two.

There does not appear to be a way to get the exact behaviour that you described, in which both tests are applied at each iteration. But it likely could be made possible if the stoppingcriterion option were to allow a list of criteria to be applied at each iteration. For example, stoppingcriterion=[function_value, relative] could allow more than one test to be applied. With that syntax then at each iteration it would first check the function value for acceptance, and only if that failed then move on to a relative error test.

The documentation could be also improved with respect to what is in effect another stopping criterion it already uses, which is maxiterations the maximal number of allowed iterations. By experiment one can show that maxiterations takes precedence over the listed stoppingcriterion options. For example,

> Student:-NumericalAnalysis:-Bisection(
>               p, x=[0,11], output=sequence,
>               stoppingcriterion=function_value, maxiterations=2); 

    [0., 11.], [0., 5.500000000], [2.750000000, 5.500000000]

So that precedence might be documented more painstakingly (read painfully?) clearly.

acer

@mattfred You wrote, "I see no reason why it should change its root finding method in any way or why it should not produce a continous smooth curve..."

 If you are going to construct the RootOf yourself, with those kind of expectations as stated and for a non-polynomial, then you might be better off using the bounding-box option rather than the root-selector option. By which I mean that you should supply some desired range, in which there is at most one root.

For example,

restart:
`ℏ`,alpha,m0,delta,`ϵg`,nu:=1,1e-24,-1e18,1,-1,2e-10:
k := (x) -> sqrt(2*m0*x*1.602/(10^19*`ℏ`^2
*(alpha+1/3*nu^2*(2/(`ϵg`+2*m0*x*1.602/(10^19*`ℏ`^2))
+1/(`ϵg`+delta+2*m0*x*1.602/(10^19*`ℏ`^2)))))):

g :=  (x, a) -> BesselJ(1, 1/1000000000*k(x)*a):

plot(g(x, 3.5), x = -5 .. 5, y = -1 .. 1);

h := (a) -> RootOf(g(x, a), x, .5):
h_alt := (a) -> RootOf(g(x, a), x, .1 .. 0.9):

Ph := plot(h(a), a = 2 .. 4, y=-2..2, color=red):
Ph_alt := plot(h_alt(a), a = 2 .. 4, y=-2..2, color=green):
plots:-display([Ph, Ph_alt]);

The basic idea is that the RootOf in h_alt(a) has a range (bounding- box) rather than a starting point (root selector). You didn't post the parameters' numeric values, so I had to fiddle them. But hopefully the code above shows (in your Maple 10- as it does in my Maple 14) that the discontinuities in the red plot can be removed as in the green plot. Of course, if you don't have a nonvarying range in which you can expect just a single root then it would get tricker, possibly with a range that depended upon `a`.

I realize that you had already figured out most of that, based upon your original Post. But I wanted to mention that supplying a range which is too narrow will often give `fsolve` difficulty. I know that sounds counter-intuitive, and it likely reflects undesirable behaviour on fsolve's part. (I believe what happens is that when fsolve's internal solving process steps "out of bounds" it immediately gives up, for that initial point. This is not good; convergence might still have attained, inside bounds, and hence that is giving up inappropriately.)

Robert mentioned alternatives to NextZero. You could construct a procedure `R` that repeatedly calls fsolve, using its 'avoid' option, to find all roots in a stated range. Then you could select the root closest to a given value. You could call `R` and supply the most recent (prior) root found. In this way you might be able to produce a continuous curve of computed roots, as the parameter varies. Here's an attempt, which seems to work,

R:=proc(f,x,rng)
   global init;
   local i, found, S, s, closest;
   S := {};
   while not(type(s,'specfunc'(anything,fsolve))) do
      s := fsolve(f,x=rng,'avoid'=S);
      if type(s,numeric) then
         S := S union {x=s};
      end if;
   end do;
   if nops(S)=0 then return 'procname'(args); end if;
   S:=sort(map(rhs,convert(S,list)));
   closest := S[1];
   for i from 1 to nops(S) do
      if abs(S[i]-init) < abs(closest-init) then
         closest := S[i];
      end if;
   end do;
   init := closest;
   closest;
end proc:

init := 0.75:
plot( 'R'(g(x,a), x, 0.1..1.0), a=2..4, y=-2..2 );

Since `plot` itself may jump around when adaptively selecting x-axis values from which to compute, it may be more prudent when using `R` to form a point-plot (possibly displayed with line style) with ordered input values than to merely call `plot` as I did above. That is due to R's reliance on the global `init` and its value as the previous computed root.

acer

There's some interesting bits to be seen in this old (2008) post on the topic.

acer

  from i from 1 to N do
     try
        res[i] := your_attempt;
     catch:
        res[i] := 0;
     end try;
  end do;

acer

Let's go back and look again at your posted code once more. Maybe there is a simple misunderstanding going on.

When you write the following, what did you intend?

     G[1] := fdiff(EIG, [1], [X[1], X[7]]);

Does `EIG` always expect two arguments, or seven arguments? If `EIG` expects seven arguments then you should call `fdiff` like this instead,

     G[1] := fdiff(EIG, [1], [X[1], X[2], X[3], X[4], X[5], X[6], X[7]]);

Note that the Vector argument `X` of `objfgradient` contains the values of the components of the 7-dimensional "point" at which you want the partial derivatives computed numerically. For `EIG` to be called correctly from within `fdiff` then all the seven numerical scalar elements of the point in question would be needed.

As an `fdiff` usage example, if one has a procedure/operator `f` that takes 5 arguments, then numerical estimation of the partial derivative with respect to the third argument -- evaluated at the given point [a,b,c,d,e] would be done with a call like fdiff(f, [3], [a,b,c,d,e]) .

I see now that my original example, of a procedure that takes only two arguments, is ambiguous. The call fdiff(f, [1], [a,b]) could just as easily be misinterpreted as only requiring the first and last argument in that list [a,b]. But in fact all entries of the given point are needed, as fdiff must pass them along when it in turn calls `f`. Sorry for any confusion.

If this is in fact the source of the problem then you'd also need to change the assignment to G[2], etc... Maybe something like this,

objfgradient := proc(X::Vector, G::Vector)
local Xlist, i, m;
   m := op(1,X);
   Xlist := convert(X,list);
   for i from 1 to m do 
      G[i] := fdiff(EIG, [i], Xlist);
   end do;
NULL;
end proc;

acer