The exponential function is entered as exp(x) not e^x in 1D Maple Notation. In that mode, `e` has no special meaning.

What you have done is solve something like the bivariate problem of maximizing b^a, but using the mere names e and x instead of b and a. And it stopped with a value of b that hits what is taken as the "infinite bound".

> with(Optimization):

> Maximize(exp(x), {x <= 1.});
Warning, no iterations performed as initial point satisfies first-order conditions

                       [2.71828182845904509, [x = 1.]]

> Maximize(exp(x), {x <= 1.}, initialpoint=[x=0.2]);

                       [2.71828182845904509, [x = 1.]]

> Maximize(e^x, {x <= 1.}); # default infinite-bound =10^20

       [                      19  [                       19        ]]
       [9.99990000000000000 10  , [e = 9.99990000000000 10  , x = 1.]]

> Maximize(e^x, {x <= 1.}, infinitebound=1e23);

       [                      22  [                       22        ]]
       [9.99989999999999982 10  , [e = 9.99990000000000 10  , x = 1.]]

You can set infolevel[Optimization]:=3 to see what value it is using for the infinite bound (and other options).

There is confusion because exp(x) gets typeset as e^x in 2D Math output (and can be also made to look like that for 2D Math input using command-completion or palette entry modes). Notice the font difference in the pretty printed output when you enter a line exactly like this (and press the enter/return key).

   exp(x), exp(1), e^x, e;

(NB. This Mapleprimes site incorrectly displays exp(x) with its "Maple Math" button as having italic "e". Bah! You have to see the distinction for exp(x) in Maple itself. You can however see exp(1) on this site, since the "Maple Math" button allows exp(1) to be displayed as as which is almost visibly distinct from which it italic like all other non-special names. It would be easier to demonstrate this here, visibly, if the display 2D math weren't rendered so awfully.)

You're hoping for something like this?

> Z:=(x,y)->g(y*exp(-x/y))*exp(x/y);

               /     /  x\\    /x\
    (x, y) -> g|y exp|- -|| exp|-|
               \     \  y//    \y/

> pdetest(u(x,y)=Z(x,y),
>         [(-x-y)*diff(u(x,y),x) - y*diff(u(x,y),y)+u(x,y)=0,u(0,y)=g(y)]);
                             [0, 0]

It looks... good. Not fantastic, but good.

If I were buying a new laptop, I'd prefer to have 4GB or more of main memory (RAM), or at least the ability to expand to that. AFAIK the VIAO VGN-FE590 series has 1GB in standard configuration, and expands to 2G only. If you are considering buying a replacement laptop then you should be able to find a VIAO with 4GB, assuming you prefer that brand. Of course, you would then want to reconsider your choice of Operating System (ie, not 32bit XP), see here.

Also, it looks like your VIAO's series has three possible models: vgnfe590a, vgnfe590b, and vgnfe590c. The a and b seem to have the Nvidia Go 7400 graphics chipset, while the c seems to have just whatever integrated chipset comes on the motherboard, ie. Intel 945GM Graphics. I don't believe that either of those configurations supports CUDA at any level. It's not really a huge deal in Maple 14, but maybe it might matter to you in the future? I think it might take 8xxxx in that graphics line, to get cuda. Just google "viao cuda" if you really want solutions here. Don't misunderstand, I don't see a reason not to to expect plotting and related Maple 14 graphics to do well with the machine you mentioned. The CUDA support in Maple 14 is restricted to a small portion of large size hardware float Matrix multiplication.

Having used the solution A1 &* A2, it is also easy to get the multiplication to be applied through by using `expand`.

note. if using `eval`, `.` might be used instead of `*` (although it doesn't matter in the scalar multiplication case).

A1:=(1/2)*I:
A2:=Matrix(2, 2, {(1, 1) = -E[0]*d[ba]/`ℏ`, (1, 2) = delta,
                  (2, 1) = -delta, (2, 2) = -E[0]*d[ab]/`ℏ`}):

T:=A1 &* A2;

                              [  E[0] d[ba]              ]
                              [- ----------     delta    ]
                       1      [    ℏ                ]
                  T := - I &* [                          ]
                       2      [                E[0] d[ab]]
                              [   -delta     - ----------]
                              [                  ℏ  ]
expand(T);

                    [  1                               ]
                    [  - I E[0] d[ba]                  ]
                    [  2                  1            ]
                    [- --------------     - I delta    ]
                    [      ℏ         2            ]
                    [                                  ]
                    [                    1             ]
                    [                    - I E[0] d[ab]]
                    [    1               2             ]
                    [   -- I delta     - --------------]
                    [    2                   ℏ    ]

eval(T,`&*`=`.`);

                    [  1                               ]
                    [  - I E[0] d[ba]                  ]
                    [  2                  1            ]
                    [- --------------     - I delta    ]
                    [      ℏ         2            ]
                    [                                  ]
                    [                    1             ]
                    [                    - I E[0] d[ab]]
                    [    1               2             ]
                    [   -- I delta     - --------------]
                    [    2                   ℏ    ]

Adding an extension to `print` can also make it look as desired. Below, there is a visual distinction between using `.` and `*`.

`print/&*` := proc(a,b) 'a . b'; end proc:

T;

                           /[  E[0] d[ba]              ]\
                           |[- ----------     delta    ]|
                   /1  \   |[    ℏ                ]|
                   |- I| . |[                          ]|
                   \2  /   |[                E[0] d[ab]]|
                           |[   -delta     - ----------]|
                           \[                  ℏ  ]/

`print/&*` := proc(a,b) 'a * b'; end proc:

T;
                          [  E[0] d[ba]              ]
                          [- ----------     delta    ]
                      1   [    ℏ                ]
                      - I [                          ]
                      2   [                E[0] d[ab]]
                          [   -delta     - ----------]
                          [                  ℏ  ]

expand(T);
                    [  1                               ]
                    [  - I E[0] d[ba]                  ]
                    [  2                  1            ]
                    [- --------------     - I delta    ]
                    [      ℏ         2            ]
                    [                                  ]
                    [                    1             ]
                    [                    - I E[0] d[ab]]
                    [    1               2             ]
                    [   -- I delta     - --------------]
                    [    2                   ℏ    ]

unassign('`print/&*`'); # undo it

T;
                            [  E[0] d[ba]              ]
                            [- ----------     delta    ]
                     1      [    ℏ                ]
                     - I &* [                          ]
                     2      [                E[0] d[ab]]
                            [   -delta     - ----------]
                            [                  ℏ  ]

if type(ans,specfunc({`+`,name=range},{sum,Sum})) then
  op(0,ans)(op([1,1],ans),op(2,ans))
  = add(op(0,ans)(-op([1,i],ans),op(2,ans)),i=2..nops(op(1,ans)));
end if;

I'm sure that there's a lot of room to make this leaner and faster.

Below, I'm not even considering whether the basic approach (which in essence boils down to the `sort` & check), is best. I also haven't considered whether the description is correct, or whether it hits the goal. I didn't see your posted code return 3/15, which you said was valid, for n=4,8,16. Did I misunderstand?

So, I'm just looking lightly at the code. I believe that these edits produce the same results as from what you posted. But you can check that.

Off the top, how about a `break` statement in the section that you called "a bit of a mess". There's no need to complete that inner loop, if you know early that it fails the requirement.

There is no need for the x[]. Just use the S[].

Inline the action that `f` does inside the code. No need for the expense of function calls. No need for the expense of calling piecewise, if the only purpose if to implement a conditional.

For example:

restart:
do_it := proc(n)
local K,S,i,reject,x,j,k;
  K := NULL:
  S := Array(1..n):
  for i from (2^n)/(2*n) + 1 to (2^n)/n - 1 by 2 do
    S[1] := i/(2^n - 1);
    for j to n-1 do
      if S[j]<=1/2 then
        S[j+1]:=2*S[j];
      else
        S[j+1] := 2*S[j]-1:
      end if;
    end do:
    S := sort(S):
    reject:=0;
    for k from 2 to n do
      if S[k] >= k/n or S[k] <= (k-1)/n then
        reject:=1:
        break;
      end if:
    end do:
    if reject=0 then
      K := K, S[1]
    end if:
  end do:
  K := [K];
end proc:

d1:d2:d3:gc():
st:=time():
do_it(16);

[2479   2539   2671   2683    959   193   1007   1009   1087   
[-----, -----, -----, -----, -----, ----, -----, -----, -----, 
[65535  65535  65535  65535  21845  4369  21845  21845  21845  

  225   3449   3557   259   261   3929   3941 ]
  ----, -----, -----, ----, ----, -----, -----]
  4369  65535  65535  4369  4369  65535  65535]

time()-st;
                             0.078

I said that I wasn't consider the actual methodology. But I'm curious, why were the two inequalities in the "check" on S[k] both strict? Did you want one strict and one nonstrict, or won't it even come up and matter?

Try inserting a line like

uses(LinearAlgebra);

right at the start of your procs that access routines from that package.

Or use the long-forn names like LinearAlgebra:-MatrixInverse throughout, for each routine from that package.

And the same goes for routines from other packages.

If you upload the full work, we could likely show the edits exactly, if you encounter problems with the fixing.

Avoid just issuing with(LinearAlgebra) outside the proc, to get it to work. I may help right now, but as you've seen it's an erratic solution (esp. if you ever intend on anyone else's using your code, or you using it a year from now when you've forgotten such subtleties: you want it solid now).

It is easy to utilize the solution returned from `solve`, by using `eval`. See examples below.

> restart:

> eqns := {x+y=a, x-y=b};
                     {x - y = b, x + y = a}

> solu := solve(eqns, {x, y});
                /    1     1          1     1  \ 
               { x = - b + - a, y = - - b + - a }
                \    2     2          2     2  / 

> eval( x/y, solu );
                           1     1   
                           - b + - a 
                           2     2   
                          -----------
                            1     1  
                          - - b + - a
                            2     2  

> eval( sin(x), solu );
                            /1     1  \
                         sin|- b + - a|
                            \2     2  /

> eqns; # still accessible, yay
                     {x - y = b, x + y = a}

> map(int,lhs~(eqns),x); # still possible, yay
                    /1  2        1  2      \ 
                   { - x  - y x, - x  + y x }
                    \2           2         / 

> neweqns := {x-y=A, x-y=B}; # still possible, yay
                     {x - y = A, x - y = B}

There are drawbacks to having x and y be assigned the values of some solution. It means that the original equations or formulae involving x and y are no longer immediately available for further symbolic manipulation. And new equations involving x and y are also not immediately attainable (you'd have make the extra effort to unassign x and y, either using the `unassign` command or by assigning them their uneval-quoted own names).

> assign(solu); # OK, now let's actually assign x and y the values

> x, y;
                     1     1      1     1  
                     - b + - a, - - b + - a
                     2     2      2     2  

> eqns; # no longer accessible, sad face
                         {a = a, b = b}

> newesteqns := {x+2*y=Q, x-2*y=R}; # not possible, sad face
                /  1     3        3     1      \ 
               { - - b + - a = Q, - b - - a = R }
                \  2     2        2     2      / 

> unassign('x'): unassign('y'): # extra work needed first

> newesteqns := {x+2*y=Q, x-2*y=R}; # possible, only after unassigning
                   {x - 2 y = R, x + 2 y = Q}

I don't understand the obsession with assigning scalar values to names. It just robs that name from the available workspace. Why not just assign equations like x=... and y=... to a list or to a name, and then utilize that with `eval`?

Partly, I blame the manuals, for encouraging such assignment to names right away. New users will latch onto anything they are shown. So they should be shown the most carefully considered answers. Using `assign` gets a slick effect, but anyone who does not already know how to utilize the results from `solve` is likely not going to know how to manage any occasionally unwanted consequences of that assignment.

You don't need to stick the f[i] into a list in order to pass them. (But of course you could do that too.)

Note that you don't have to place a type-check on the parameters in the procedure definition. See the first example below.

In consequence, there are many ways to do this. A sample are below.

> restart:
> # with no prior assignment of f as Vector, list, etc,
> # then these next lines make f a table.
> f[1] := x[2]+3*x[4]+x[5]-x[6]=2:
> f[2] := x[1]+2*x[2]-x[3]+x[4]+3*x[6]+x[7]=13:
> f[3] := x[1]+x[2]-3*x[3]-x[4]+x[5]-x[8]=26:

> p:=proc(F)
>   op(F);
> end proc:

> p([f[1]]);
                x[2] + 3 x[4] + x[5] - x[6] = 2

> p([f[1],f[2],f[3]]);
      x[2] + 3 x[4] + x[5] - x[6] = 2, 
      x[1] + 2 x[2] - x[3] + x[4] + 3 x[6] + x[7] = 13, 
      x[1] + x[2] - 3 x[3] - x[4] + x[5] - x[8] = 26

> p:=proc(F::list(equation))
>   op(F);
> end proc:

> p([f[1]]);
                x[2] + 3 x[4] + x[5] - x[6] = 2

> p([f[1],f[2],f[3]]);
      x[2] + 3 x[4] + x[5] - x[6] = 2, 
      x[1] + 2 x[2] - x[3] + x[4] + 3 x[6] + x[7] = 13, 
      x[1] + x[2] - 3 x[3] - x[4] + x[5] - x[8] = 26

> p:=proc(F::table)
>   # F[2] is available here, or we can get them all...
>   entries(F,'nolist');
> end proc:

> p(f);
      x[2] + 3 x[4] + x[5] - x[6] = 2, 
      x[1] + 2 x[2] - x[3] + x[4] + 3 x[6] + x[7] = 13, 
      x[1] + x[2] - 3 x[3] - x[4] + x[5] - x[8] = 26

> p:=proc(F::seq(equation))
>   F;
> end proc:

> p(f[1]);
                x[2] + 3 x[4] + x[5] - x[6] = 2

> p(f[1],f[2],f[3]);
      x[2] + 3 x[4] + x[5] - x[6] = 2, 
      x[1] + 2 x[2] - x[3] + x[4] + 3 x[6] + x[7] = 13, 
      x[1] + x[2] - 3 x[3] - x[4] + x[5] - x[8] = 26

Here is an attempt

> pol2cart:=proc(theta,rho)
>    if type(theta,scalar) and type(rho,scalar) then
>       rho*cos(theta),rho*sin(theta),`if`(nargs=3,args[3],NULL);
>    else
>       zip(`*`,rho,map(cos,theta)), zip(`*`,rho,map(sin,theta)),
>       `if`(nargs=3,args[3],NULL);
>    end if;
> end proc:

> pol2cart(Pi/6,5);

                                 5  (1/2)  5
                                 - 3     , -
                                 2         2

> pol2cart(<Pi/6,Pi/4,Pi/2>, <2,3,4>);

                           [  (1/2) ]  [   1    ]
                           [ 3      ]  [        ]
                           [        ]  [3  (1/2)]
                           [3  (1/2)]  [- 2     ]
                           [- 2     ], [2       ]
                           [2       ]  [        ]
                           [        ]  [   4    ]
                           [   0    ]            

> pol2cart(Pi/6,5,10);

                               5  (1/2)  5    
                               - 3     , -, 10
                               2         2    

> pol2cart(<Pi/6,Pi/4,Pi/2>, <2,3,4>, <7,8,9>);
                         [  (1/2) ]  [   1    ]     
                         [ 3      ]  [        ]  [7]
                         [        ]  [3  (1/2)]  [ ]
                         [3  (1/2)]  [- 2     ]  [8]
                         [- 2     ], [2       ], [ ]
                         [2       ]  [        ]  [9]
                         [        ]  [   4    ]     
                         [   0    ]                 

You can now get the following behaviour:

> eval(parse(Matlab:-FromMatlab("pol2cart(THETA,RHO)",string)));
                       RHO cos(THETA), RHO sin(THETA)

But I wasn't able to get the full expression sequence returned, without that 'string' option. This is unfortunately only half the desired result.

> Matlab:-FromMatlab("pol2cart(THETA,RHO)");
Evaluating:
pol2cart(THETA, RHO);

                               RHO cos(THETA)

I also tried that last example, after using Matlab:-AddTranslator but it didn't make a difference.

I didn't try using it with any FromMFile example.

The code above deliberately avoids (~) elementwise operators, in order to work in versions before Maple 13. For those who like a laugh...

      if parse(StringTools:-Split(
                  kernelopts(version)," ")[2][1..-2]) >= 13 then
         eval(parse(`_Rho *~ (cos~(_Theta)), _Rho *~ (sin~(_Theta))`),
              [_Rho=rho,_Theta=theta]);
      else
         zip(`*`,rho,map(cos,theta)), zip(`*`,rho,map(sin,theta));
      end if;

I've noticed that the FromMatlab and FromMFile environment is a bit more of a Matlab simulator for Maple than a translator. Is that good, or bad? It looks like herclau is saying that it's easier to translate source entirely by hand than it is to take incomplete and broken results from FromMatlab or FromMFile and repair them. Have I got that close to right? Is this because the Matlab "translator" produces more high-level routine calls, instead of lower level translation? Is it weaker than FromMma?

Regardless of whether sets or tables might be better for whatever task the OP has, here's a comparison.

restart:
X:=[$1..10^5]:
t:=[$1..10^3]:
st,ba,bu:=time(),kernelopts(bytesalloc),kernelopts(bytesused):
eval(subs(`=`~(t,nULL),X),nULL=NULL):
time()-st,kernelopts(bytesalloc)-ba,kernelopts(bytesused)-bu;
                    0.344, 1441528, 1249108

restart:
X:=[$1..10^5]:
t:=[$1..10^3]:
st,ba,bu:=time(),kernelopts(bytesalloc),kernelopts(bytesused):
remove(z->type(z,identical(op(t))),X):
time()-st,kernelopts(bytesalloc)-ba,kernelopts(bytesused)-bu;
                   2.578, 10156220, 10067364

restart:
X:=[$1..10^5]:
t:=[$1..10^3]:
st,ba,bu:=time(),kernelopts(bytesalloc),kernelopts(bytesused):
subsindets(X,identical(op(t)),z->NULL):
time()-st,kernelopts(bytesalloc)-ba,kernelopts(bytesused)-bu;
                   4.954, 48225664, 798240352

restart:
X:=[$1..10^5]:
t:=[$1..10^3]:
st,ba,bu:=time(),kernelopts(bytesalloc),kernelopts(bytesused):
Y:=X:
for i from 1 to nops(t) do Y:=subs(t[i]=NULL,Y) end do:
time()-st,kernelopts(bytesalloc)-ba,kernelopts(bytesused)-bu;
                   5.765, 49536144, 796094064

At much smaller sizes (of both X and t) they all get quicker. But the first above (Preben's Method 1) is still fastest of those shown.

A list is probably the wrong thing to be using, if you're going to be doing a lot of such augments.

As Preben mentioned, the issue for lists is that augmenting them entails an (inefficient) copying of the object. The same is true of merely changing a single entry of a list (while keeping the total number of entries the same).

Preben showed in his response that you could use a table instead. But you can see that it's a little unwieldy, as far as the syntax goes.

In recent versions of Maple, round-bracket () indexing allows augmentation of Vectors, Matrices, and Arrays. You can "grow" a Vector with this simple syntax. And I believe it can be shown by experiment that Maple might be augmenting with slightly more space than is requested (for non-hardware datatypes) so that a few repeated augments can be done with only a single initial copy occuring up front.

> c := <1,2>;
                                       [1]
                                  c := [ ]
                                       [2]

> c(3) := Pi;
                                       [1 ]
                                       [  ]
                                  c := [2 ]
                                       [  ]
                                       [Pi]

> c := Vector[row]([1,2]);
                                 c := [1, 2]

> c(3) := exp(1);
                             c := [1, 2, exp(1)]

And, of course, you can alter a pre-existing entry of a Vector without any copying occuring since a Vector is a mutable data structure.

No, rand() produces uniformly distributed samples and your posted Statistics package example involves normally distributed samples, so their distributions are very different.

Of course, Statistics also allows for uniformly distributed samples as well. So you could compare that against rand().

But note that rand() is just using RandomTools, as its help page describes:

  rand calls RandomTools[MersenneTwister][GenerateInteger] or RandomTools[MersenneTwister][NewGenerator]
  depending on whether or not a number or procedure is to be returned.  It is more efficient
  to make these calls directly than to call rand.

•  The random number generator used by rand can be seeded by using the randomize or
  RandomTools[MersenneTwister][SetState] functions.

  The algorithm used by rand in Maple versions up to and including 9.5 has been moved
  into the RandomTools package as RandomTools[LinearCongruence].

The biggest difference is efficiency. If you need big samples from a continuously distributed random variable then Statistics:-Sample is the leaner way to go. It produces double precision hardware floats (in datatype=float[8] rtables) by default.

Have you read John's series of posts on this topic? (especially here, here). He talks about and compares those methods mentioned in the ?rand help page.

What I don't see anywhere obvious is a description of what Statistics uses underneath (for Uniform samples, say). Erik had lots to say about how custom distributions are handled. But I didn't see anything stating whether uniform samples produced by the Statistics package also use whatever method those produced by the RandomTools package use (or not).

You could try alpha=0, beta=1, to resolve two of those parts, which may lead you to sigma=1/2 to get continuity. But then the second deriv tends to -1/3 as x->0 from the right. Too bad.

You could try to match the second deriv at x=0 and as x->0 from the right. If that limit must be zero, then that may lead you to s=0 or s=1/6. But neither of those will allow matching from the left (if alpha and beta must be not merely real but also constant).

I'm sure that someone will let us know, if I messed this up.

restart:
# First way: Deduce alpha=0, and then deduce sigma to match
# second deriv. Then test that against 0th and 1st derivs.
ee:=(1+a*x^2)^ln(x);
                                  ln(x)
                        /       2\     
                        \1 + a x /     
eval(ee,a=0);
                               1
ff:=ln(1+s*x^2)/(1-cos(x));
                            /       2\
                          ln\1 + s x /
                          ------------
                           1 - cos(x) 
limit(diff(ff,x,x),x=0);
                           1        2
                           - s - 2 s 
                           3         
solve(%=0);
                                 1
                              0, -
                                 6
limit(eval(ff,s=0),x=0); # darn, we want 1 here
                               0
limit(eval(ff,s=1/6),x=0); # darn, we want 1 here
                               1
                               -
                               3
# Second way: find sigma to match function value (0th deriv)
# and then test against 1st and 2ns derivs.
limit(ee,x=0); # true also with alpha=0
                               1
solve(limit(ff,x=0)=%, s);
                               1
                               -
                               2
FF:=eval(ff,s=%);
                            /    1  2\
                          ln|1 + - x |
                            \    2   /
                          ------------
                           1 - cos(x) 
limit(diff(FF,x),x=0);
                               0
limit(diff(FF,x,x),x=0); # darn, we want 0 here
                               -1
                               --
                               3 

@icedragon

Another issue that we didn't mention before: you have y=-20..20 so you better use that instead of y(x) which is another beast entirely.

What you have now is not even an equation. That is, there's no left-hand-side = right-hand-side.

What you have instead is just this:

  (x^((1/(3)))* (y(x)^(2)+ 5* y(x)) + 3 *(x^(2)+y(x)^(2))*sin(xy(x)),  i^()

What you might want is something more like this (an equation):

  (x^((1/(3)))* (y(x)^(2)+ 5* y(x)) + 3 *(x^(2)+y(x)^(2))*sin(x*y(x)) = i^(3)

Try this below. (Exactly as it appears. Cut and paste. Without a space between with and (plots) too, nor after implicitplot and its bracketed arguments.).

restart:
with(plots):

col := ["Red", "Green", "Blue", "Orange", "Teal",
        "Azure", "BlueViolet", "Yellow", "DarkOliveGreen", "DarkRed"]:

for i from 1 to 10 do
 z[i]:=implicitplot(x^((1/(3)))* (y^(2)+ 5* y)
                    + 3 *(x^(2)+y^(2))*sin(x*y) = i^(3),
                    x = -10...10,  y= -20...20,
                    color = col[i], numpoints = 2000)
end do:
display(z[3]);
display(seq(z[i], i = 1 .. 10));

animate(implicitplot,
        [x^((1/(3)))* (y^(2)+ 5* y)
         + 3 *(x^(2)+y^(2))*sin(x*y) = 'i'^(3),
         x = -10...10,  y= -20...20, numpoints = 2000],
        'i'=1..10);

The last, red, plot that this should show above is an animation. If you single-click to select it then you should see Play/FF/Reverse buttons in the menubar.