What's the objection to using piecewise?

> f := proc(x) piecewise(x>=0 and x<1,1,0); end proc:
> int(f(x),x=0..2);
                                       1
 
> int(f,0..2);
                                       1
 
> f := piecewise(x>=0 and x<1,1,0):
> int(f,x=0..2);
                                       1

Consider your original,

> f:=proc(x) if x>=0 and x<1 then 1 else 0 fi; end proc:
> f(x);
Error, (in f) cannot determine if this expression is true or false: 0 <= x and
x < 1

Note that int(f(x),...) or plot(f(x),...) will try to evaluate f(x) right away. That's why you got an error. Because at that moment, x is not assigned and so the comparison produces that error. You could also use the operator form of calling int (or plot), even with your original definition of f.

> int(f,0..2);
                                       1

acer

Actually, you set up an equation Determinant(PL) = 0, and started working with that. So, either deal with lhs() of that equation, or don't set it up as an equation.

EqPl:=simplify(Determinant(PL)=0):
coeffs(lhs(EqPl), [x, y, z]);

altern := simplify(Determinant(PL)):
coeffs(altern, [x, y, z]);

BTW, that subs() call didn't assign to anything. Did you intend it to do so?

acer

The sort of membership that you are considering involves checking whether your candidate Matrix is identical to any member of the set A.

But Matrices are mutable data structures (their entries can change). Two mutable structures are not considered identical just because they happen to have the same entries at one given moment. That's because changing the entry of one of them wouldn't affect the other.

The same happens for direct comparison of Matrices. Consider,

> evalb( Matrix(2,2,[[0,0],[0,0]]) = Matrix(2,2,[[0,0],[0,0]]) );
                                     false
 
> X := Matrix(2,2,[[0,0],[0,0]]):
> Y := Matrix(2,2,[[0,0],[0,0]]): # not same object as X
> evalb( X = Y );
                                     false
 
> X[1,1]:=17:
> Y[1,1];
                                       0
 
> Y := X: # same object as X
> evalb( X = Y );
                                      true

> X[1,1]:=13:
> Y[1,1];
                                       13

For your set A, you could instead compare entries something like this,

> ormap(t->EqualEntries(t,Matrix(2,2,[[0,0],[0,0]])),A);
                                      true
 
> ormap(t->EqualEntries(t,Matrix(2,2,[[1,1],[1,1]])),A);
                                      true

You can also use more sophisticated comparisons. Eg, here is a floating point comparison (allowing for  only 1 ulp difference  when compared using 3 digits)

> ormap(t->verify(t,Matrix(2,2,[[1,1],[1,0.999]]),
>                 'Matrix'('float'(1,digits=3))),A);
                                      true
 
> ormap(t->verify(t,Matrix(2,2,[[1,1],[1,0.998]]),
>                 'Matrix'('float'(1,digits=3))),A);
                                     false

Having introduced verify, one could do your example like this,

> verify(Matrix([[1,1],[1,1]]),A,'member'('Matrix'));
                                     true
 
> verify(Matrix([[1,1],[2,2]]),A,'member'('Matrix'));
                                     false

acer

> L := 0.5*m*diff(x(t),t)^2 - m*g*x(t):

> diff(frontend(diff,[L,diff(x(t),t)]),t) = frontend(diff,[L,x(t)]);
                                  / 2      \
                                  |d       |
                            1.0 m |--- x(t)| = -m g
                                  |  2     |
                                  \dt      /

acer

> T:=convert(x^n*exp(x),FormalPowerSeries,x);
                                 infinity
                                  -----    (k + n)
                                   \      x
                            T :=    )     --------
                                   /         k!
                                  -----
                                  k = 0

 > Z:=value(subs(infinity=5,T));
                              (2 + n)    (3 + n)    (4 + n)    (5 + n)
              n    (1 + n)   x          x          x          x
        Z := x  + x        + -------- + -------- + -------- + --------
                                2          6          24        120
 
> S:=series(eval(Z,n=3),x,8);
                     3    4        5        6         7      8
               S := x  + x  + 1/2 x  + 1/6 x  + 1/24 x  + O(x )
 
> value(subs(n=3,T));
                                    3
                                   x  exp(x)
> series(%,x,8);
                   3    4        5        6         7      8
                  x  + x  + 1/2 x  + 1/6 x  + 1/24 x  + O(x )

> Z:=value(subs(infinity=20,T)):
> MultiSeries:-series(Z,x,5);
                          (2 + n)    (3 + n)    (4 + n)
          n    (1 + n)   x          x          x             (5 + n)
         x  + x        + -------- + -------- + -------- + O(x       )
                            2          6          24

acer

I suspect that this may explain how djc obtained this form,

   evalm((`&*`(`&*`(A, B), 2))*B-`&*`(B, Id));

If the ?evalm help-page, being viewed, has its Examples toggled to 2D Math (using the small button symbol at right-end of iconic menubar), then pasting the relevant line to a Worksheet as 1D Maple input produces the following,

   evalm((`&*`(`&*`(A, B), 2))*B-`&*`(B, Id))

But it the Examples are toggled to 1D Maple input, then it gets copied and pasted into a Worksheet as the following 1D Maple input,

   evalm((A &* B) &* (2*B)-B &* Id);

I also tried pasting as 2D Math, with both levels of Typesetting (Extended and Standard). But who can tell, what is really lying "underneath" displayed 2D Math? I don't like the idea that 2D Math might be WYSINWYH (what you see is not what you have).

acer

By default, units are not combined upon addition.

> restart:
> 3*Unit(microliters)+4*Unit(liters);
                                3 [uL] + 4 [L]

They are combined, if the Units:-Standard package is loaded. The default volume unit for results is meter^3.

> restart:
> with(Units:-Standard):
> 3*Unit(microliters)+4*Unit(liters);
                                 4000003   [ 3]
                                ---------- [m ]
                                1000000000

One can change it, so that microliters is the new default for volume in results.

> restart:
> with(Units:-Standard):
> Units:-UseSystem(Units:-AddSystem(SI,uL));
> 3*Unit(microliters)+4*Unit(liters);
                                 4000003 [uL]

acer

> T := table([ (1,1)=a, (2,2)=s ]);
                     T := table([(1, 1) = a, (2, 2) = s])

> Matrix(2,2,{op(op(eval(T)))});
                                   [a    0]
                                   [      ]
                                   [0    s]

> B[1,2]:=y:
> B[2,1]:=z:
> Matrix(3,3,{op(op(eval(B)))});
                                 [0    y    0]
                                 [           ]
                                 [z    0    0]
                                 [           ]
                                 [0    0    0]

acer

It is the principal third root. There are three cube roots of -1, lying equally spaced around the unit circle in the complex plane.

> evalc( (-1)^(1/3) );
                                            1/2
                               1/2 + 1/2 I 3
 
> convert(%,polar);
                                          Pi
                                polar(1, ----)
                                          3

For a non-principal root,

> surd(-1,3);
                                      -1

There is a short, related, comment in the sqrt help-page (4th bullet item in Description). Compare with exp(1/3 * ln(-1)).

acer

It's not enough to simply assign to the entries B[i,j]. You also have to create B as a matrix or Matrix. Otherwise, B will get created as a table, which is a different structure and not acceptable to matrixplot.

> restart:

> B[1,1]:=4:
> whattype(B);
                                    symbol
 
> type(B,table);
                                     true
 
> plots:-matrixplot(B);
Error, (in plots/matrixplot) invalid input: convert/Matrix
expects its 1st argument, M, to be of type
{Array, Matrix, Vector, array, list}, but received B

> B:=matrix(1,1,[]):
> B[1,1]:=4:
> whattype(B);
                                    symbol
 
> type(B,array);
                                     true
 
> plots:-matrixplot(B); # it works

acer

> f := x -> sin(x):
> RootFinding:-NextZero(f,0.0);
                                  3.141592654
 
> df := x -> evalf(D(f)(x)):
> RootFinding:-NextZero(df,0.0);
                                  1.570796327

> signum( evalf(D[1,1](f)(%)) );
                                      -1

So, the next zero of df, moving positively from 0.0, is at 1.57... and it's a local maximum of f since the second derivative of f is negative at that point. That might be the fastest way when coded efficiently, if just the single next extremum is wanted. If the higher derivative test was not showing the correct sign for what you wanted (maximum or minimum) then NextZero could be re-applied to df starting from the current point, and the test re-checked. Repeat until an extremum of the sort requested is found. It should be pretty straightforward to write a procedure which does all this.

If the function application, f(x), is a univariate polynomial with all other parameters instantiated at known numeric values then fsolve(diff(f(x),x)) would produce all the real roots of the first derivative at once. (Or Student:-Calculus1:-ExtremePoints might be used.) Those could then be sorted, and the higher derivative test could be used to show whether each is maximum or minimum, or to find the smallest positive extremum of the requested type.

acer

It looks as if the solver is getting trapped at a saddle point, depending on the initial value from which it begins solving. (I'm not sure if it matters, that the point lies on a constraint or intersection of constraints. Perhaps the first feasible point tried is that unfortunate point, although if that's the case then it's a pity something else isn't subsequently tried.)

One workaround might be to symbolically simplify the constraints so that it picks a different initial point which allows convergence to some other point . Another workaround may be to provide a judicious (and feasible) initial point directly to the solver.

> x[max] := 10:
> x[min] := 0:
> y[max] := 10:
> y[min] := 0:
> delta := 1:
> b:=2:
> d:=5:
> expr:=((2*c+a-d-2*b)*delta-(1/2)*c+(1/2)*d-b+a)^2/(9*b-9*a);
                                  /3 c             \2
                                  |--- + 2 a - 17/2|
                                  \ 2              /
                          expr := -------------------
                                       18 - 9 a
 
> cons:={b-a >= c-d,
>        (2*b-a-2*c+d)*delta <= 2*b-2*a-2*c+2*d,
>         c-d-b+a <= (2*b-a-2*c+d)*delta};
    cons := {c <= 7 - a, -a - 2 c <= 5 - 2 a - 2 c, c + a <= 16 - a - 2 c}
 
> bounds:= a = x[min] .. b, c = d .. y[max];
                       bounds := a = 0 .. 2, c = 5 .. 10
 
> solA := Optimization:-NLPSolve(expr,cons,bounds,maximize);
Warning, no iterations performed as initial point satisfies
 first-order conditions
               solA := [-0., [a = 0.500000000000000000, c = 5.]]
 
> eval(diff(expr,a),op(2,solA));
                                      0.
 
> eval(diff(expr,c),op(2,solA));
                                      0.
 
> solve(cons union {a>=0, a<=2, c>=5, c<=10});
                                                         2 a
                  {0 <= a, 5 <= c, a <= 1/2, c <= 16/3 - ---}
                                                          3
 
> Optimization:-NLPSolve(
>           expr,
>           solve(cons union {a>=0, a<=2, c>=5, c<=10}),
>           bounds,maximize);
                   [0.0555555555555555525, [a = 0., c = 5.]]

> Optimization:-NLPSolve(expr,cons,bounds,maximize,
>                        initialpoint=[a=0.4,c=5]);
                   [0.0555555555555555525, [a = 0., c = 5.]]

As far as plotting goes, you can get nice visuals using the interactive assistant (toggling the "maximize" checkbox, then Solve button, then Plot button). Compare the plots in the assistant for these two calls,

Optimization:-interactive(expr,cons,a=0..1/2,c=5..16/3);
Optimization:-interactive(expr,cons,a=0..1,c=5..16/3);

acer

What is this shorthand for? Does it stand for a procedure call, or a procedure definition?

   w := proc (Q1, Q2) -->evalf function....

Hopefully, it's not a procedure being (re)created inside a double-loop.

If you post all the code, there's a better chance someone might optimize it well.

acer

Are you saying that, for either subset, you just want the its own members to be linearly independent? If so, then try to find just two pairs. Ask yourself, what does it mean for a pair of vectors to be linear independent? How can you tell quickly that a pair of vectors are linearly dependent?

Or, are you requiring something about both sets at once (a direct sum of the two subspaces spanned by the members of the two sets, say)? What would it mean if a 5x4 matrix, with four of the five vectors as its columns, had a rank of 4?

acer

> restart:
> W := Matrix([[1,0],[0.1,1]]);
                                     [ 1     0]
                                W := [        ]
                                     [0.1    1]
 
> map(limit, LinearAlgebra:-MatrixFunction(W,x^t,x), t=infinity);
                            [      1.           0.]
                            [                     ]
                            [Float(infinity)    1.]

You might also consider looking up the spectral radius of a Matrix. (Especially the first theorem on that page.)

 acer