Are you using a full colon, after the outermost `end do`, to suppress output?

acer

> restart:

> e:= G*(.7500000000*B*sqrt((1.-2.*a+a^2)*B^2/(a^2+B^2-2.*B^2*a+B^2*a^2))-.7500000000*B*sqrt((1.-2.*a+a^2)*B^2/(a^2+B^2-2.*B^2*a+B^2*a^2))*a+.7500000000*a*sqrt(a^2/(a^2+B^2-2.*B^2*a+B^2*a^2))) + (1-G)*(-6.250000000*10^(-12)*sqrt(a^2/(a^2+B^2-2.*B^2*a+B^2*a^2))*(5.833333333*10^9*B*sqrt((1.-2.*a+a^2)*B^2/(a^2+B^2-2.*B^2*a+B^2*a^2))-7.500000000*10^9)*(4.*a-3.*B*sqrt((1.-2.*a+a^2)*B^2/(a^2+B^2-2.*B^2*a+B^2*a^2))+3.*B*sqrt((1.-2.*a+a^2)*B^2/(a^2+B^2-2.*B^2*a+B^2*a^2))*a-3.*a*sqrt(a^2/(a^2+B^2-2.*B^2*a+B^2*a^2)))^2-6.250000000*10^(-12)*(5.833333333*10^9*sqrt(a^2/(a^2+B^2-2.*B^2*a+B^2*a^2))-7.500000000*10^9)*B*sqrt((1.-2.*a+a^2)*B^2/(a^2+B^2-2.*B^2*a+B^2*a^2))*(4.-4.*a-3.*B*sqrt((1.-2.*a+a^2)*B^2/(a^2+B^2-2.*B^2*a+B^2*a^2))+3.*B*sqrt((1.-2.*a+a^2)*B^2/(a^2+B^2-2.*B^2*a+B^2*a^2))*a-3.*a*sqrt(a^2/(a^2+B^2-2.*B^2*a+B^2*a^2)))^2+0.3645833333e-1*sqrt(a^2/(a^2+B^2-2.*B^2*a+B^2*a^2))*B*sqrt((1.-2.*a+a^2)*B^2/(a^2+B^2-2.*B^2*a+B^2*a^2))*(4.-3.*B*sqrt((1.-2.*a+a^2)*B^2/(a^2+B^2-2.*B^2*a+B^2*a^2))+3.*B*sqrt((1.-2.*a+a^2)*B^2/(a^2+B^2-2.*B^2*a+B^2*a^2))*a-3.*a*sqrt(a^2/(a^2+B^2-2.*B^2*a+B^2*a^2)))^2+5.625000000*10^(-11)*(1.000000000*10^10+5.833333333*10^9*sqrt(a^2/(a^2+B^2-2.*B^2*a+B^2*a^2))*B*sqrt((1.-2.*a+a^2)*B^2/(a^2+B^2-2.*B^2*a+B^2*a^2))-7.500000000*10^9*B*sqrt((1.-2.*a+a^2)*B^2/(a^2+B^2-2.*B^2*a+B^2*a^2))-7.500000000*10^9*sqrt(a^2/(a^2+B^2-2.*B^2*a+B^2*a^2)))*(-1.*B*sqrt((1.-2.*a+a^2)*B^2/(a^2+B^2-2.*B^2*a+B^2*a^2))+B*sqrt((1.-2.*a+a^2)*B^2/(a^2+B^2-2.*B^2*a+B^2*a^2))*a-1.*a*sqrt(a^2/(a^2+B^2-2.*B^2*a+B^2*a^2)))^2):

> # first, get rid of those suspiciously simple looking floats.
> # (Isn't there as easier way to do this, btw??)

> f:=simplify(subsindets(e,float,t->10^10*identify(10^(-20)*(identify(identify(t)*10^10))))) 
>    assuming a>0, a<1, B>0, B<1, G>0, G<1:

> f:=simplify(f,size);
                          /  /                              
            1             |  |3         // 4   2   2    \  4
------------------------- |3 |- (G - 1) ||B  - -- B  + 1| a 
                    (3/2) \  \4         \\     27       /   
  /       2  2    2\                                        
4 \(a - 1)  B  + a /                                        

     /    4   4   2\  3   /   4   2   2\  2      4      4\ 
   + |-4 B  + -- B | a  + |6 B  - -- B | a  - 4 B  a + B | 
     \        27   /      \       27   /                 / 

                    (1/2)                                         
  /       2  2    2\        //     2\  2      2      2\ //        
  \(a - 1)  B  + a /      + \\1 + B / a  - 2 B  a + B / \\(G - 1) 

                                                             \\
   2        \  3   /           2    \  2    2              2\||
  B  - G + 1/ a  + \(3 - 2 G) B  + G/ a  + B  (G - 3) a + B /||
                                                             //

> ## That attains its (global w.r.t B and a) maximum when a=0, for any given G.
> ## If in doubt, I suppose one could plot it.
> ## It also attains a global max (w.r.t B and a) for all B along a=1, when G=1. See below.
> #plot3d(eval([e],G=1.0),a=0..1,B=0..1,axes=boxed); # or any G>0, G<1

> f_at_a_equals_zero := simplify(limit(f,a=0),size) assuming G>=0, G<=1, B>=0, B<=1;
                     1             2   3  
                     -- (9 G - 9) B  + - B
                     16                4  

> maximize(f_at_a_equals_zero,B=0..1,G=0..1,location); # global max for all B and G, a=0
                   3   /[                3]\ 
                   -, { [{B = 1, G = 1}, -] }
                   4   \[                4]/ 

But the goal was to find the global maximum for each value of G, not the maximum over all G.

At G=1 the global maximum (in terms of a and B) also occurs all along a=1, and in fact it is constant there,

> eval(f,[a=1,G=1]);
                               3
                               -
                               4

> optB := solve(diff(f_at_a_equals_zero,B),B);
                                2    
                          - ---------
                            3 (G - 1)

> simplify(eval(f_at_a_equals_zero,B=optB)); # the answer
                                 1    
                          - ---------
                            4 (G - 1)

> optG_func := unapply(%,G):

> optG_func(0.15);
                          0.2941176470

> optG_func(0.10);
                          0.2777777778

> optG_func(1/3),optG_func(1.0/3);
                        3              
                        -, 0.3750000000
                        8              

acer

This answer below relates to your follow-up challenge which seems quite different fromn your original post.

Please let us know if you think of other special conditions, as yet unmentioned.

It seems that you are trying to `select` based on the presence of the word "how" and not just the substring "how" so that "thisishowIlikeit" would not match. So, you could Split the strings at delimiter " ", or you could Search for " how " and check both ends as well.

> aa:=[["he",45,123,76,1.0,"4"],["know how",4,9,34,"3.2","5"]]:

> select(t->ormap(z->member("how",StringTools:-Split(z," ")),indets(t,string)),aa);
              [["know how", 4, 9, 34, "3.2", "5"]]

NB. You mentioned StringTools:-Has, before. But that is likely not what you want, since it is a character locater not a substring matcher.

NB. Regarding your earlier talk about parse, the idea of parsing strings containing spaces, to rely on 2D implicit multiplication, looks dubious at best. What about "this 1is 2how $I like it, eh" or "some are 1-1 and some are many-1".

acer

> simplify((combine(expr1-expr2,power))) assuming nonnegative;
                               0

There were some spam responses (now flagged). Is the whole thread just some link-bait spam?

acer

The problem of a missing space (hence, no actual implicit multiplication) also occurs for `λl`.

That same 2D Math mistake persists in this other post (a duplicate, IMO).

acer

The issue is that the extracted sub-Matrix has global n, and not the assumed, local n.

You can fix it up in several ways. You can create a whole new Matrix, a replica of the extrcated sub-Matrix but with entries evaluated so as to have the assumed n (using copy, or the Matrix constructor, or eval). Or you can use the 'inplace' option to rtable_eval and evaluate the entries of the extracted sub-Matrix without creating yet another Matrix.

Note that merely extracting the sub-Matrix causes creation of a new Matrix. You are thus trying to avoid creating yet another Matrix.

The result of atataa() already has the assumed, local n. There is no need, or effect gained, from copying or evaluating it.

> restart:

> M:=Matrix([[n]]):
> N:=Matrix([[n]]):

> atataa:=proc(a) local s;
>   s:=Matrix(1,1);
>   s[1,1]:=a[1,1];
>  s;
> end proc:

> with(LinearAlgebra):

> assume(n,positive);

> Equal(atataa(M), SubMatrix(N,1..-1,1));
                             false

> Equal(atataa(M), copy(eval(SubMatrix(N,1..-1,1)))); # another Matrix created
                              true

> Equal(atataa(M), Matrix(N[1..-1,1])); # another Matrix created
                              true

> Equal(atataa(M), N[1..-1,1]); # alternate syntax for getting the sub-Matrix
                              false

> subm_N:=SubMatrix(N,1..-1,1): # or N[1..-1,1], the same deal
> Equal(atataa(M), subm_N);
                             false

> rtable_eval(subm_N,inplace): # no Matrix created
> Equal(atataa(M), subm_N);
                              true

A note on efficiency. You wrote that this Equal test will occur inside some loops. Consider altering atataa() to accept both V and Vtmp as a pair of Matrix arguments. You could then re-use the very same Vtmp each time. Even if you have to pay the cost of zeroing Vtmp at the start of atataa (using ArrayTools:-Fill, say) it might still be more efficient than producing garbage Matrix objects inside atataa (which must be collected via memory-management) with each iteration.

acer

I got an answer in two different ways, using the Optimization package. But I did have to make a few corrections (similar to what was done here for the DirectSearch approach).

One thing that I did not do (but which the code shouts for, I think) is improvements to make the objective evaluate more quickly. Even a single evaluation of the objective at a given numeric point can be very slow. And optimization methods can request a lot of objective evalutations. For example, there are some numeric integrals, embedded deep in the mess (sorry) of the objective procedure, which are identically repeated several times. It's likely taking at least five times longer than it should have to, and that's not even considering whether there is some purely exact way to get a single expression for any of its current fsolve & integrals. I also didn't test to see whether the Digits:=30 was really needed, or best, for the fsolve inside procedure SX.

Because each objective evalutation is so slow, and given that you know the objective is smooth (and possibly even monotonic in each variable), then I would expect that a local optimization approach (even if performed by some global-capable optimization package) could be faster than a full strength global optimization method.

But first, a note. That error message about no iteration steps take due to first order conditions already being met, from Optimization, is often due to an internal bug related to a problem with computing numeric derivatives. It can be side-stepped by either using the unconstrained Nelder-Mead nonlinearsimplex method (slow convergence rate here, and many objective evaluations) with which I got a result. Or it can be side-stepped by supplying the obective in so-called "Matrix form", along with a hand-built routine to compute its numerical derivatives.

# I used this redefinition of SX, because sometimes the fsolve call fails.
SX := proc (egx, egy)::set; local j, sx, SEqn; global Eg1, Eg2;
 SEqn := {seq(evalf(subs(subs(Eg1 = egx, Eg2 = egy, EV), Eqn(i))), i = 1 .. n)};
 Digits := 30; sx := fsolve(SEqn, X);
 if type(sx, specfunc(anything, fsolve)) then
   return {seq(x[i] = 1, i = 1 .. n)};
 else return sx; end if;
end proc

# And I used this redefinition of PEffSQ
PEffSQ := proc (Egx, Egy) local res;
  if Egy <= Egx then
    res := 100*evalf(subs(subs(Eg1 = Egx, Eg2 = Egy, EV), SX(Egx, Egy), PCellMax/PInc));
  end if;
  if `not`(type(res, realcons)) then return 0; else return res; end if;
end proc

I noted that the Egy<=Egx inside PEffSQ could be turned into a constraint, but didn't go that route even though the above changes make the objective nondifferentiable at the boundary.

At that point, I could do the following workaround

> objf := proc(V::Vector)
>   PEffSQ(V[1],V[2]);
> end proc:

> objfgradient := proc(X::Vector,G::Vector)
>   G[1] := fdiff( PEffSQ, [1], [X[1],X[2]] );
>   G[2] := fdiff( PEffSQ, [2], [X[1],X[2]] );
>   NULL;
> end proc:

> Optimization:-NLPSolve( 2, objf, 'objectivegradient'=objfgradient,
>         'initialpoint'=Vector([2.3, 1.4]), 'method'='sqp', 'maximize' );

                  [                     [ 2.24196645821080  ]
                  [48.3648680038370032, [                   ]
                  [                     [ 1.41837902813424 ]]

I have a hope that any problem in NLPSolve could be fixed, so that it could automatically use the above or some other way to compute numerical derivatives and succeed more often for "Operator" form arguments.

acer

The evalf/Int approach should use sophisticated techniques. It also allows use to add an accuracy tolerance as an extra epsilon=value option inside that inert Int() call `igral`. In contrast, the Student:-NumericalAnalysis methods are more for teaching principals: stock (almost "toy") methods without optional accuracy options.

You could simply raise the working precision and recompute, and then see which of the two earlier results was a better approximation to the more accurate high-Digits result. The simplest way to do that is to make an assignment such as Digits:=20 or so.

When you call evalf(Int(...)) as in your first approach with igral, Maple automatically adjusts the target accuracy (epsilon) based on Digits the working precision ebvironment variable. So when you raise Digits and recompute, the result should become more accurate. Here, that's just simpler than adding the epsilon=value option into the Int() call.

And there is another way to get the same effect, of raising the working precision. That is to change the evalf to evalf[d] in your first approach. That makes Maple use `d` for Digits for the given computational argument to evalf.

Also, you can set infolevel[`evalf/int`]:=2 (note that was all lowercase) and see more information about what evalf/Int is doing and what accuracy tolerance it is using, etc.

acer

The call to MatrixInverse must somehow bring along a connection to the name LinearAlgebra, because otherwise how would Maple know to put that in the error string?

The Original Poster's pair of commands obvously work OK when done alone in a fresh session, so the challenge becomes to find a sequence of commands which could be done neforehand, so that the given pair of lines would produce the error.

First let's consider a sequence of commands that works as I'd expect,

> restart:
> with(Student):

> with(LinearAlgebra):
> MatrixInverse( Matrix([[4]]) );
                                            MatrixInverse([4])

That makes sense. When I loaded the Student package, it rebound the name LinearAlgebra since that is one of the Student package's exports. So the subsequent command with(LinearAlgebra) was actually loading Student:-LinearAlgebra. And there is no export `MatrixInverse` in the Student:-LinearAlgebra subpackage, so the name MatrixInverse didn't match any known procedure, and what gets returned is an unevaluated function call. So far, so good.

Now another example, whose final two lines are closer to but not exactly the same as what was originally reported,

> restart:
> with(Student):

> with(:-LinearAlgebra):
> MatrixInverse(Matrix([[4]]));
Error, MatrixInverse is not a command in the LinearAlgebra package

Now I am surprised. I expected with(:-LinearAlgebra) to load by force the top level LinearAlgebra package and not any subpackage export from the Student package. Maybe it did. But, if so, then why did the MatrixInverse call fail? It's as if Maple figured out that the MatrixInverse call should correctly try to invoke :-LinearAlgebra:-MatrixInverse but the "LinearAlgebra" bit of some lookup then found the "wrong" name on account of the rebinding to Student:-LinearAlgeba. I am a little surprised, and maybe it's a bug.

Now another example, whose last two lines exactly match the original example,

> restart:
> with(LinearAlgebra):
> with(Student):

> with(LinearAlgebra):
> MatrixInverse(Matrix([[4]]));
Error, MatrixInverse is not a command in the LinearAlgebra package

I suspect that this too might be unintentional and have the same cause as the middle example, and hence be a bug. It does seem undesirable behaviour, and if these are not bugs then how can they be clearly documented (especially now that the GUI allows "easy" package loading)?

acer

By using eval, one can substitute the prior solutions involving x and y (as a set of equations) without having to actually assign to x and y. This is nice because it allows one to create new systems involving x and y as unassigned variable names, even after solving the initial system for x and y.

> restart:

> eqxy1:=x+2*y=4;
                          x + 2 y = 4

> eqxy2:=x^2+y^2=16;
                           2    2     
                          x  + y  = 16

> xysol:=solve({eqxy1,eqxy2},Explicit);
                               /    -12      16\ 
              {x = 4, y = 0}, { x = ---, y = -- }
                               \     5       5 / 

> xysol[1];
                         {x = 4, y = 0}

> xysol[2];
                       /    -12      16\ 
                      { x = ---, y = -- }
                       \     5       5 / 

> eqrsxy1:=s*x+y^2*r=11;
                               2       
                        s x + y  r = 11

> eqrsxy2:=r^2+s^2=x*y;
                          2    2      
                         r  + s  = x y

> eval({eqrsxy1,eqrsxy2},xysol[1]);
                    /           2    2    \ 
                   { 4 s = 11, r  + s  = 0 }
                    \                     / 

> eval({eqrsxy1,eqrsxy2},xysol[2]);
             /  12     256          2    2   -192\ 
            { - -- s + --- r = 11, r  + s  = ---- }
             \  5      25                     25 / 

> solve(eval({eqrsxy1,eqrsxy2},xysol[1]),Explicit);
           /    11        11\    /      11        11\ 
          { r = -- I, s = -- }, { r = - -- I, s = -- }
           \    4         4 /    \      4         4 / 

> solve(eval({eqrsxy1,eqrsxy2},xysol[2]),Explicit);
           /    4400     3             (1/2)  
          { r = ---- + ----- I 15164737     , 
           \    4321   17284                  

                  4125     16             (1/2)\    /
            s = - ----- + ----- I 15164737      }, { 
                  17284   21605                /    \

                4400     3             (1/2)  
            r = ---- - ----- I 15164737     , 
                4321   17284                  

                  4125     16             (1/2)\ 
            s = - ----- - ----- I 15164737      }
                  17284   21605                / 

acer

The fact that Maple appears to be doing the exact integral is a bit of an illusion, because all those individual erf() calls in the expanded "result" would likely have to be approximated numerically at some point.

You can consider doing numeric integration directly, instead of using value() to get an apparent exact result. Consider,

> restart:

> expr:=(a-b*delta)*(c-delta)*
>   1/(sigma*sqrt(2*Pi))*exp(-(delta-mu)^2/2/sigma^2):

> myvals:=[a=20,b=2,c=1,sigma=99,mu=1/10]:

> igral:=Int(expr,delta=0..a/b):

> evalf(subs(myvals,igral)); # better
                         -0.9386112295

> evalf(eval(value(igral),myvals)); # worse
                         -0.9386105310

>evalf[20](eval(value(igral),myvals)); # check
                    -0.93861122946668528680

The "better" answer above did pure numeric integration on the original expression. The second, "worse" answer resulted from applying value() and getting an intermediary answer with lots of erf() calls, and only subsequently applying evalf() to get a float answer.

Of course, one can raise Digits (or call evalf[d](...) to get the same effect) and attain greater accuracy. And that's what the final "check" did. But it's also worthwhile noticing which of the two initial attempts was more accurate (better) for this example.

There are some cases where the ability to transform a single integral into a (set of) calls to some named special function is a good thing (because direct numeric integration might be otherwise hard or notoriously innaccurate when compared to a specialized internal routine designed to evaluate the special function numerically). But as the above example illustrates, that's not always the case.

acer

The help-page ?ExpectedValue could be more clear about this usage, especially because ExpectedValue(X,u) does not appear to work properly for X a RandomVariable instead of a discrete Sample.

restart:
with(Statistics):

Q:=RandomVariable(Distribution(PDF=(t->f(t)))):

ExpectedValue(Q);

ExpectedValue(U(Q));

Q:=RandomVariable(Normal(mu,sigma)):

ExpectedValue(U(Q));

And one could replace U above, as desired. There's no need for creation of a new, weighted random variable proper.

For example, this is just a wrong result (and looks like a bug to me, since it's more in line with the documented calling sequences than is the last call above),

> restart:
> with(Statistics):
> Q:=RandomVariable(Normal(mu,sigma)):
> ExpectedValue(Q,t->1-exp(-phi*t));
                                       mu

acer

There are two problems in the initial 2D Math input of e1.

The factor that looks like  is actually (a-`b&delta;`) instead of (a-b*delta). I suggest retyping it in full (and don't just try and insert a space or *).

Also it is missing either a space or an explicit multiplication symbol between the first two factors in brackets. That makes it a function application rather than a multiplication of two bracketed subexpressions. I suggest retyping it with an explicit *.

So, it has instead of either or .

Once changed, Maple can do the integral directly, without expansion by hand, using lowercase int() instead of evalf(Int(...)). That produces the result containg erf() calls.

If you intend on doing numeric integration with evalf(Int(..)) then you'll need to use numeric values for the unknowns a,b,c,sigma, and mu in advance of the call.

If you encounter poor accuracy from with numeric integration or with floating-point evaluation of the exact symbolic integration result (containing erf calls) then you could try increasing the environment variable Digits. See ?Digits.

Alternatively, is this something you could more simply compute using the Statistics package?

If you continue to have difficulties with 2D Math input then you could consider switching from a Document to a Worksheet and from 2D Math input to 1D Maple input.

acer

The error message occurs because the procedure has a break statement in it (2nd line from the end) which is not inside any for..do...end loop. Try removing it.

What are you trying to accomplish with that break? Did you intend that the procedure would instead exit (return) at the line with that invalid break? If so, then you could use an error statement instead of merely printf'ing that error message. Ie, replace,

  printf("Error: At least 2 planes do not have cut point. Please change the plane equation, and try again.");
  break;

with this,

  error "At least 2 planes do not have cut point. Please change the plane equation, and try again.";

acer

> res:=evalc(int(1/x,x));
                            /1   1          \   
                ln(|x|) + I |- - - signum(x)| Pi
                            \2   2          /   

> convert(res,piecewise);
                        {  ln(-x) + Pi I          x < 0
                        {
                        { ln(x) + 1/2 I Pi        x = 0
                        {
                        {      ln(x)              0 < x

> convert(res,abs);
                               /1     x  \   
                   ln(|x|) + I |- - -----| Pi
                               \2   2 |x|/   

I believe that the value for the condition x=0 in the above piecewise may be affected by the environment settings for signum(0) and Heaviside(0).

Also discussed here and here, including G A Edgar's suggestion,

> Re(int(1/x,x));
                          ln(abs(x))

From the help-page ?ln

      For complex-valued expressions x, ln(x) = ln(abs(x))+I*argument(x), where
       -Pi<argument(x)<=Pi.   Throughout Maple, this computation is taken to be the
      definition of the principal branch of the logarithm. 

acer