The default order of elements in a Matrix is Fortran_order, in which matrix elements listed as

A:=Matrix(2,symbol=a);

                           [a[1, 1]    a[1, 2]]
                      A := [                  ]
                           [a[2, 1]    a[2, 2]]

seq(A(i),i=1..4);

                  a[1, 1], a[2, 1], a[1, 2], a[2, 2]

seq(i,i=A);

                  a[1, 1], a[2, 1], a[1, 2], a[2, 2]

In this order, to get what you need, the matrix has to be transposed,

A^%T;

                         [a[1, 1]    a[2, 1]]
                         [                  ]
                         [a[1, 2]    a[2, 2]]

If you defined the C_order in the matrix, the transposing wouldn't be needed,

A:=Matrix(2,symbol=a, order=C_order);

                           [a[1, 1]    a[1, 2]]
                      A := [                  ]
                           [a[2, 1]    a[2, 2]]

seq(A(i),i=1..4);

                  a[1, 1], a[1, 2], a[2, 1], a[2, 2]

seq(i,i=A);

                  a[1, 1], a[1, 2], a[2, 1], a[2, 2]

Alec

" For Maple the best way to sum over an array would be to use add."

I wouldn't be that sure about that. Anyway, even using add, it can be done faster.

n := 10^8:
A := LinearAlgebra:-RandomVector( n, outputoptions=[datatype=integer[8]] ):
time[real](add( A[i], i=1..n )); 
                                         
                                29.718

time[real](add( i, i=A ));

                                19.297

Converting A to `+` is about as fast as that,

time[real](convert(A,`+`));

                                19.282

The following doesn't work in this example,

time[real](`+`(entries(A,nolist)));
Error, object too large

but for smaller n it works, and it is faster,

n := 10^7:
B := LinearAlgebra:-RandomVector( n, outputoptions=[datatype=integer[8]] ):
time[real](add( B[i], i=1..n )); 

                                2.933

time[real](convert(B,`+`));

                                1.919

time[real](`+`(entries(B,nolist)));

                                0.779

The following is even faster,

time[real](rtable_scanblock(B,[],Sum));

                                0.483

and it works in the original example as well,

time[real](rtable_scanblock(A,[],Sum));

                                4.820

The following also doesn't work in the original example,

time[real](trunc(Statistics:-Mean(A)*10^8));
Error, (in Statistics:-Mean) not enough memory to allocate rtable

but in the smaller example where it works, it is faster,

time[real](trunc(Statistics:-Mean(B)*10^7));

                                0.124

Alec

It is also possible to write the printf commands as

printf("%{}15a\n%15f",data[1,1],data[2]);
              t            v(t)            x(t)
       0.000000        0.000000        1.000000
       0.100000       -0.099833        0.995004

Alec

It is also possible to write the printf commands as

printf("%{}15a\n%15f",data[1,1],data[2]);
              t            v(t)            x(t)
       0.000000        0.000000        1.000000
       0.100000       -0.099833        0.995004

Alec

MakeKeyMatrix is very inefficient, it needs about a second to produce a 100 by 100 matrix

tt:=time(): to 100 do A:=MakeKeyMatrix(100) od: 
time()-tt;

                                98.671

Here is about 1000 times faster procedure,

g:=proc(n) 
    local A,r; 
    uses LinearAlgebra:-Modular; 
    r:=0; 
    while r=0 do 
        A:=Create(127,n,n,'random',('integer')[]);
        r:=Determinant(127,A) od; 
    A 
end: 

tt:=time(): to 100 do A:=g(100) od: time()-tt;

                                0.109

Alec

MakeKeyMatrix is very inefficient, it needs about a second to produce a 100 by 100 matrix

tt:=time(): to 100 do A:=MakeKeyMatrix(100) od: 
time()-tt;

                                98.671

Here is about 1000 times faster procedure,

g:=proc(n) 
    local A,r; 
    uses LinearAlgebra:-Modular; 
    r:=0; 
    while r=0 do 
        A:=Create(127,n,n,'random',('integer')[]);
        r:=Determinant(127,A) od; 
    A 
end: 

tt:=time(): to 100 do A:=g(100) od: time()-tt;

                                0.109

Alec

By the way, there is a concavity option in Student:-Calculus1:-FunctionChart (or FunctionPlot) drawing arrows up on the intervals where a function is stritly convex (i.e. concave up) and arrows down where a function is strictly concave.

Alec

Acer,

Thank you! I extremely value all yours, Robert Israel's, Joe Riel's, Jacques Carette, Dave Linder, and several other people posts.

If Maplesoft hired me 5 or 7 years ago, Maple would be already much better by now, I think.

Alec

I agree and it is great that you started doing that. The procedure as it is should work if the is could be corrected.

That's an old topic - what 1/x - 1/x should return. Currently it is 0 in Maple, but that creates problems as one can see in the last example.

That is one thing where domains may be useful - if 1/x is an element of the field of rational functions, it should be 0. If it is an expression representing a function, it should be undefined at 0.

Alec

You could post it either here, or in the Mathematica newsgroup. There are a lot of skilled people there as well.

Alec

You could post it either here, or in the Mathematica newsgroup. There are a lot of skilled people there as well.

Alec

Linear functions are convex, just not strictly convex. Same as constants are increasing, but not strictly increasing.

Alec

Fails for abs though,

isconvex(abs(x));
isconvex:   falling back to the definition

                                 FAIL

isconvex(abs(x^3));
isconvex:   falling back to the definition

                                 FAIL

And while replacing abs(x) with sqrt(x^2) gives the correct answer,

isconvex(sqrt(x^2));

                                 true

the following is wrong,

isconvex(-sqrt(x^2));

                                 true

That was caused by the following Maple feature

diff(sqrt(x^2),x,x);

                              2
                             x          1
                         - ------- + -------
                             2 3/2     2 1/2
                           (x )      (x )

normal(%);

                                  0

is(%%<=0);
                                 true

is(%%%=0);

                                 true

Alec

Yes, and works also in Maple 12 Classic. But it was undocumented at that time. Actually, even now it is not mentioned in the ?comment help page in Classic Maple; it is mentioned in ?comment in Standard though.

Alec

Yes, and works also in Maple 12 Classic. But it was undocumented at that time. Actually, even now it is not mentioned in the ?comment help page in Classic Maple; it is mentioned in ?comment in Standard though.

Alec