acer

32622 Reputation

29 Badges

20 years, 44 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

As shown in response to your duplicate post, you can specify a global solver for Optimization:-Maximize.

In this example, it may happen that the default local solver used by Maximize also returns the global maximum. But you can also force the method, with an option.

> with(Optimization):

> nx := Maximize(abs(.5547114632*x^2+1.130478381*x+.9883691714-exp(x)),
               x = -1 .. 1, method=branchandbound);

            nx := [0.0456887233851495, [x = -0.436328905794407]]

acer

For some reason, Maximize is using its "local" solver on this problem by default.

> restart:

> with(Optimization):

> infolevel[Optimization]:=1:

> nx := Maximize(
>         abs(.5539395590*x^2+1.130864333*x+.9891410756-exp(x)),
>         x=-1.0..1.0); # "local solver"

NLPSolve: calling NLP solver
NLPSolve: using method=quadratic
NLPSolve: number of problem variables
            nx := [0.0452333936742979, [x = -0.438621263028586]]

> nx := Maximize(
>         abs(.5539395590*x^2+1.130864333*x+.9891410756-exp(x)),
>         x=-1.0..1.0,method=branchandbound); # "global solver"

NLPSolve: calling NLP solver
NLPSolve: using method=branchandbound
NLPSolve: number of problem variables
             nx := [0.0454683311488107, [x = 0.560938672728197]]

acer

The command from the Optimization package is Maximize, not maximize.

> with(Optimization):
> nx := Maximize(abs(.5547114632*x^2+1.130478381*x+.9883691714-exp(x)),
>                x = -1 .. 1);

            nx := [0.0456887233851494, [x = -0.436328906042811]]

> restart:
> nx := Optimization:-Maximize(
>                abs(.5547114632*x^2+1.130478381*x+.9883691714-exp(x)),
>                x = -1 .. 1);

            nx := [0.0456887233851494, [x = -0.436328906042811]]

acer

Capitalization matters in Maple. The command from the Optimization package is Maximize, not maximize.

> with(Optimization):

> nx := Maximize(abs((29.83428688*x-57.16914857)/(x-2.196093020)-exp(x)),
>                x = 5.0 .. 6.0);
                     nx := [371.399468551146, [x = 6.]]

> restart:
> Optimization:-Maximize(abs((29.83428688*x-57.16914857)/(x-2.196093020)-exp(x)),
>                x = 5.0 .. 6.0);
                        [371.399468551146, [x = 6.]]

acer

Since you've now seen the `expand` command, here's some more Maple that might help you,

> L := [-1/2,2/3,7/3,4+3*I,4-3*I];

[-1 2 7 ]
[--, -, -, 4 + 3 I, 4 - 3 I]
[2 3 3 ]

> facp := mul((x-r), r in L);

/ 1\ / 2\ / 7\
|x + -| |x - -| |x - -| (x + (-4 - 3 I)) (x + (-4 + 3 I))
\ 2/ \ 3/ \ 3/

> p:= expand(facp);

5 21 4 175 811 3 373 2 29
x - -- x + --- + --- x - --- x - -- x
2 9 18 6 6

[solve(p)];

[-1 7 2 ]
[--, -, -, 4 - 3 I, 4 + 3 I]
[2 3 3 ]

and `fsolve` is an alternative to `solve` that you could find useful for higher degree univariate polynomials.

acer

It can get trickier, if there are more "mismatches" at the range ends. A procedure can be constructed which uses a single intermediate that seems to do the trick (in at least some more cases).

> combminus:=proc(s::specfunc(anything,Sum),t::specfunc(anything,Sum))
>    local temp;
>    if op(1,s) = op(1,t) then
>       temp:=Sum(op(1,s),i=op(1,rhs(op(2,s)))..op(2,rhs(op(2,t))));
>       value(combine(s-temp)+combine(temp-t));
>    else
>       s-t;
>    end if;
> end proc:


> s1 := Sum(x(i), i = n-1 .. m-1):
> s3 := Sum(x(i), i = n+1 .. m+1):

> value(combine(s3-s1));

              /  m - 1         \   /  m + 1      \
              | -----          |   | -----       |
              |  \             |   |  \          |
              |   )            |   |   )         |
              |  /      (-x(i))| + |  /      x(i)|
              | -----          |   | -----       |
              \i = n - 1       /   \i = n + 1    /

> combminus(s3,s1);

               x(m) + x(m + 1) - x(n - 1) - x(n)

# The single intermediate s2 (as in DuncanA's Answer) does not suffice alone.
# (Not that anyone has claimed it would, of course.)
> s2 := Sum(x(i), i = n .. m):

> value(combine(s3-s2)
>       +combine(s2-s1));

     /  m + 1      \   /  m          \   /  m - 1         \
     | -----       |   |-----        |   | -----          |
     |  \          |   | \           |   |  \             |
     |   )         |   |  )          |   |   )            |
     |  /      x(i)| + | /    (-x(i))| + |  /      (-x(i))|
     | -----       |   |-----        |   | -----          |
     \i = n + 1    /   \i = n        /   \i = n - 1       /

          /  m       \
          |-----     |
          | \        |
          |  )       |
        + | /    x(i)|
          |-----     |
          \i = n     /

# No need to introduce this many intermediates,
# as `combminus` solved it above.
> s2_a := Sum(x(i), i = n .. m-1):
> s2_b := Sum(x(i), i = n+1 .. m-1):
> s2_c := Sum(x(i), i = n+1 .. m):

> value(combine(s3-s2_c)
>       +combine(s2_c-s2_b)
>       +combine(s2_b-s2_a)
>       +combine(s2_a-s1));

               x(m) + x(m + 1) - x(n - 1) - x(n)

# Another example.
> s1 := Sum(x(i), i = n-1 .. m+1):
> s3 := Sum(x(i), i = n .. m-1):

> value(combine(s3-s1));

               /  m + 1         \   /m - 1     \
               | -----          |   |-----     |
               |  \             |   | \        |
               |   )            |   |  )       |
               |  /      (-x(i))| + | /    x(i)|
               | -----          |   |-----     |
               \i = n - 1       /   \i = n     /

> combminus(s3,s1);

                  -x(m) - x(m + 1) - x(n - 1)

acer

with(Student:-NumericalAnalysis):
de:=diff(y(t),t)=-exp(3*y(t)),y(0)=1:
P1:=RungeKutta(de,t=1,submethod=rk4,output=plot,numsteps=10,
               plotoptions=[color=green]):
P2:=Euler(de,t=1,output=plot,numsteps=10,
          plotoptions=[color=black]):
plots:-display([P1,P2]);

Maybe someone can show how to make the legend better, showing the colored symbols.

acer

The problem is that the units in numerator and denominator do not simplify and cancel by default.

> a/(b+t*Unit(kg));

                      7 Unit('kg')            
             -----------------------------
             10 Unit('kg') + t Unit('kg')

You can make them "cancel" in several ways.

> a:=7*Unit(kg):
> b:=10*Unit(kg):

> y:=t->ln(a/(b+t)):
> plot(combine(y(t*Unit(kg)),units));

> y:=proc(t) uses Units:-Standard;
>            ln(a/(b+t)); end proc:
> plot(y(t*Unit(kg)));

> y:=t->ln(a/(b+t)):
> plot(simplify(y(t*Unit(kg))));

acer

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

First 285 286 287 288 289 290 291 Last Page 287 of 339