Interesting question.

Smart money says that someone could improve the following,

> PE := A.X+B.U+L.(Y-C.X-D.U);
            A . X + B . U + L . (Y - C . X - D . U)

> rx:=proc(x,z)
  local ox,sx,i;
  if type(x,`*`) then
    ox,sx:=selectremove(t->type(t,`.`),[op(x)]);
    ox:=op(ox);
  else
    ox,sx:=x,1;
  end if;
  if type(ox,`.`) and op(nops(ox),ox)=z then
      [`*`(op(sx),`.`(seq(op(i,ox),i=1..nops(ox)-1))),z];
  else x;
  end if;
end proc:

> rightcollect:=proc(x,z)
local spl,t,locx;
  locx:=subs(`&*`=`.`,expand(subs(`.`=`&*`,x)));
  if type(x,`+`) then
    spl:=selectremove(t->type(t,list),map(rx,[op(locx)],z));
    `+`(op(spl[2]))+`+`(seq(op(1,t),t in spl[1])).z;
  else x;
  end if;
end proc:

> rightcollect(PE,X);
          B . U + L . Y - L . D . U + (A - L . C) . X

> PredictorEstimator := diff(X(t), t) = A.X(t)+B.U+L.(Y-C.X(t)-D.U);
   d                                                          
  --- X(t) = A . (X(t)) + B . U + L . (Y - C . (X(t)) - D . U)
   dt                                                         

> map(rightcollect,PredictorEstimator,X(t));
   d                                                         
  --- X(t) = B . U + L . Y - L . D . U + (A - L . C) . (X(t))
   dt                                                        

A related item: is it just an oversight that `expand` can deal with `&*` but not with `.`? Does anyone else have trouble writing an `expand/.` routine to do x->subs(`&*`=`.`,expand(subs(`.`=`&*`,x)))? Is it anything to do with `expandon`?

acer

inner_proc:=proc(M::integer,
                 A::Vector(datatype=float[8]),
                 B::Vector(datatype=float[8]))
local count::integer, i::integer;
count:=0;
for i from 1 to M do
  if sqrt(A[i]^2+B[i]^2) < 0.5 then
    count:=count+1;
  end if;
end do;
return 4.0*count/M;
end proc:

cinner_proc:=Compiler:-Compile(inner_proc):

p:=proc(N) # N number of points
local F,X,Y,n,inner_proc,P;
  n:=N;
  F:=Statistics:-RandomVariable(Uniform(-1/2, 1/2));
  X:=Statistics:-Sample(F,n);
  Y:=Statistics:-Sample(F,n);
  cinner_proc(n,X,Y);
end proc:

st:=time(): p(2*10^5); time()-st;
                      3.14218000000000020
                             0.063

I tried to resist the strong temptation to change how it works, and to be faithful to the given method.

acer

inner_proc:=proc(M::integer,
                 A::Vector(datatype=float[8]),
                 B::Vector(datatype=float[8]))
local count::integer, i::integer;
count:=0.0;
for i from 1 to M do
  if sqrt(A[i]^2+B[i]^2) < 0.5 then
    count:=count+1.0;
  end if;
end do;
return 4.0*count/M;
end proc:

cinner_proc:=Compiler:-Compile(inner_proc):

p:=proc(N) # N number of points
local F,X,Y,n,inner_proc,P;
  n:=N;
  F:=Statistics:-RandomVariable(Uniform(-1/2, 1/2));
  X:=Statistics:-Sample(F,n);
  Y:=Statistics:-Sample(F,n);
  cinner_proc(n,X,Y);
end proc:

st:=time(): p(2*10^5); time()-st;
                      3.14250000000000007
                             0.078

acer

What you are doing now is a bit like this (simplified to just 1 loop),

> restart:

> lc:=1:
> for ii from 1 to 3 do
>   T[ii]:=i->lc;
>   lc:=lc+1;
> end do:

> T[1](), T[2](), T[3]();
                                    4, 4, 4
 
> lc:=17:
> T[1](), T[2](), T[3]();
                                  17, 17, 17

But it seems like you want behaviour more like this,

> restart:

> lc:=1:
> for ii from 1 to 3 do
>   T[ii]:=unapply(lc,i);
>   lc:=lc+1;
> end do:

> T[1](), T[2](), T[3]();
                                    1, 2, 3
 
> lc:=17:
> T[1](), T[2](), T[3]();
                                    1, 2, 3

Is that closer to what you are after?

Notice also that there is no need to have the T indices be the same as the formal parameters of the operators (ie, i,j,k). Having it so may be just additionally confusing to you.

acer

Sorry, but 100*cos(20) is an exact quantity. (Recall, Maple is a computer algebra system, in which there is an important distinction since algebraic manipulation of exact quanitites is useful.) You can apply evalf or evalhf to it, to get a floating-point quantity.

Or, you could ensure that the argument to sin is itself a float, so that it returns a float. This is probably a better solution for you than explicit using explicit evalf/evalhf calls, if you intend to use Compile on the procedure.

Maple will not automatically cast an exact numeric quantity to a float when assigning into a float[8] Matrix, except in the cases that the exact quantity is an integer or a rational.

acer

> with(LinearAlgebra):
>
> A:=Matrix([[A11,A12],[A21,A22]]):
>
> InvA := Matrix([[(A11-A12&*A22^(-1)&*A21)^(-1),
>           -A11^(-1)&*A12&*(A22-A21&*A11^(-1)&*A12)^(-1)],
>         [-A22^(-1)&*A21&*(A11-A12&*A22^(-1)&*A21)^(-1),
>          (A22-A21&*A11^(-1)&*A12)^(-1)]]);
InvA :=
 
    [             1
    [--------------------------- ,
    [      //        1 \       \
    [A11 - ||A12 &* ---| &* A21|
    [      \\       A22/       /
 
     // 1        \                 1             \]
    -||--- &* A12| &* ---------------------------|]
     |\A11       /          //        1 \       \|]
     |                A22 - ||A21 &* ---| &* A12||]
     \                      \\       A11/       //]
 
    [ // 1        \                 1             \
    [-||--- &* A21| &* ---------------------------| ,
    [ |\A22       /          //        1 \       \|
    [ |                A11 - ||A12 &* ---| &* A21||
    [ \                      \\       A22/       //
 
                 1             ]
    ---------------------------]
          //        1 \       \]
    A22 - ||A21 &* ---| &* A12|]
          \\       A11/       /]
 
>
> map(normal, A . subs(`&*`=`*`,InvA) );
                                   [1    0]
                                   [      ]
                                   [0    1]
 
>
> data:={A11=RandomMatrix(2,2), A12=RandomMatrix(2,2),
>        A21=RandomMatrix(2,2), A22=RandomMatrix(2,2)}:
>
> A:=Matrix(eval([[A11,A12],[A21,A22]],data)):
>
> InvA:=Matrix(eval(subs(`&*`=`.`,
>           [[(A11-A12&*A22^(-1)&*A21)^(-1),
>             -A11^(-1)&*A12&*(A22-A21&*A11^(-1)&*A12)^(-1)],
>            [-A22^(-1)&*A21&*(A11-A12&*A22^(-1)&*A21)^(-1),
>             (A22-A21&*A11^(-1)&*A12)^(-1)]] ),data) ):
>
> A . InvA;
                              [1    0    0    0]
                              [                ]
                              [0    1    0    0]
                              [                ]
                              [0    0    1    0]
                              [                ]
                              [0    0    0    1]

acer

@Hamidreza Is there a timing difference if you compute

add(add(evalf(f[i, j](i, j, x, y)), i = 1 .. 50), j = 1 .. 10);

instead of

evalf(add(add(f[i, j](i, j, x, y), i = 1 .. 50), j = 1 .. 10));

Also, could you use evalhf instead of evalf, in the former of those two?

And what is D1()? Is it always the same, or does it change for each i and j (using globals)? If unchanging, then can you pull it out of the summation as a multiplicative constant, or evalf(D()) up front? That too could make a difference.

And, taking another tack, how about first finding Q:=sum(...,i=1..N) to get a formula, finding Q50:=eval(Q,N=50), and then computing add(evalf(Q50),j=1..100) or similar?

If you raise Digits, and then recompute, you can get some evidence that those relatively small imaginary components are merely artecfacts of doing the computation with floating-point numbers (ie. as an approximation).

For this example, you might be able to convince yourself of this with an assignment of only Digits:=16.

acer

Your procedure `inertia` expects to find the global xyzc to be an indexable object with three scalar values. That's how  `inertia` uses xyzc. And it works fine, the first time you call it.

But then you assign to xyzc (at the top level) something which is a list of three equations. This occurs when you compute Equate(cog(...),[xc,yc,zc]), the result of which you assign to global xyzc. So then xyzc is no longer in a form for which some of the commands internal to `inertia` are valid or permissible. The three entries of xyzc are now equations rather than scalar values. So the next invocation to `inertia` fails with an error, going wrong when it tries to square the entries of xyzc.

For example,

> (4.99=xc)^2;
Error, (in simpl/reloprod) invalid terms in product

acer

Must a local optimization (least squares) method find the global optimum for any nonlinear fit?

You could try a few things, such as setting it up yourself as a least squares problem, and then using different starting points or using a global optimizer.

I'd like to see Maple automatically detect the presence of GlobalOptimization, and use it within NonLinearFit if found. (Using yet another global optimizer might be possible but could get tricky, as it would have to accept the same syntax as Optimization and also somehow be registered for detection.)

acer

Thanks Alec! That's better for Primes, but worse for me, as the inability-to-post problem looks to be at my end.

A related question: is there anyone for whom the 2DMath displayed in Alec's response looks good? For me, it doesn't even look as good as the 2DMath in the Online Help using the exact same browser and session. It's small, and more jagged than even lack of antialising usually gives.

acer

Give it the "bug" tag, Alec.

It's an outright bug that assumptions are put on x merely by calling `series`, present in 13.01 (and probably 13.02 too). Which makes it present in the version most people are using at this date.

"Du musst Amboss oder Hammer sein." - Goethe

acer

Does LinearAlgebra:-SylvesterSolve do this?

You started with this,

> Z.C.X + C = D;

                              (Z . C . X) + C = D

So, now if you right-multiply by the inverse of X, you get,

> Z.C.X.X^(-1) + C.X^(`-1`) = D.X^(`-1`);

                                        -1           -1

                       (Z . C) + (C . (X  )) = D . (X  )

Now, if you read the help-page for LinearAlgebra:-SylvesterSolve, it says that it solves Matrix equations of the form A . X + isgn * X . B = scale*C.

And you seem to have that form, except that your C is their X, your Z is their A, your X^(-1) is their B, and your D.X^(-1) is their C. If that's correct, then how about calling,

LinearAlgebra:-SylvesterSolve(Z, X^(-1), D.X^(-1) );

You'd want to check that I didn't make a mistake in all that. If it doesn't work, maybe you could post or upload your actual Matrices.

acer

Who thinks that the original student's numeric task cannot be implemented in Maple within a factor of four or better of the Excel timing? How about a factor of 2?

Nick, please upload the Maple program if you are so permitted. If it cannot be achieved, then at least people will know more about what needs improving.

acer

Good question. When you pass an rtable (ie. any of Matrix, Array, or Vector) as an argument to a procedure call it gets passed similarly to what you called "by reference".

The procedure gets the same rtable object, not a copy, of what was passed. Hence changes made to the entries of that rtable parameter, while inside that procedure, also change the entries in the object at the higher level.

acer