There are issues with 1) automatically detecting oscillatory integrands and consequent automatic determination of the appropriate method, and 2) automatic determination of the number of subintervals into which to split the range.

Mma used to fail on some similar problems, and issue a message about the user's needing to specify a larger number of subintervals. Getting rid of that need may be what has improved in Mma 8.

In Maple 14, one can do it this way (similar in princpal to Mma's decribed older behaviour),

> st:=time():
> evalf(Int(cos(log(x)/x), x = .0001 .. 1));
                          0.3247340117
> time()-st;
                             2.094

> evalf(Int(cos(log(x)/x), x = .0001 .. 1, method=_d01akc));
Error, (in evalf/int) NE_QUAD_MAX_SUBDIV:
  The maximum number of subdivisions has
  been reached: max_num_subint = 500

> st:=time():
> evalf(Int(cos(log(x)/x), x = .0001 .. 1, maxintervals=2000, method=_d01akc));
                          0.3247340117
> time()-st;
                             0.156

The two things to see above is that I had to specify the method which I knew was "best" for the stated problem. And I had to manually increase the number of subintervals. Note that it was not always possible to even raise the subinterval number in Maple. See this old thread.

If Wolfram has figured out a reliable way to do it automatically, quickly, then... well done. There's thus room for further improvement in Maple too.

There's also room for improvement in Maple in the area of oscillatory numeric integrals over semi-infinite regions, which is a form that can also occur with a change of variables....

> Y:=Int(cos(log(x)/x), x = 0 .. 1);

                         /1              
                        |      /ln(x)\   
                        |   cos|-----| dx
                        |      \  x  /   
                       /0                

> student[changevar](y=1/x,Y,y);

                    /infinity    /  /1\  \   
                   |          cos|ln|-| y|   
                   |             \  \y/  /   
                   |          ------------ dy
                   |                2        
                  /1               y         

> op(1,%);

                             /  /1\  \
                          cos|ln|-| y|
                             \  \y/  /
                          ------------
                                2     
                               y      

> plot(%,y=1..30);

Look at how that plot decays. Imagine a scheme for oscillatory integrals for semi-infinite ranges which could quickly handle the above. Who knows, maybe that's what Mma 8 is actually doing?

acer

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